From 22ca3897e88c2cb7d23161d1f48a41c4db29abd8 Mon Sep 17 00:00:00 2001 From: Sergio Pascual Date: Wed, 14 Feb 2018 18:19:19 +0100 Subject: [PATCH] Convert methods in base scientific recipe to functions --- megaradrp/processing/extractobj.py | 256 +++++++++++++++ megaradrp/recipes/calibration/lcbstdstar.py | 30 +- megaradrp/recipes/calibration/mosstdstar.py | 10 +- megaradrp/recipes/scientific/base.py | 338 +------------------- 4 files changed, 293 insertions(+), 341 deletions(-) create mode 100644 megaradrp/processing/extractobj.py diff --git a/megaradrp/processing/extractobj.py b/megaradrp/processing/extractobj.py new file mode 100644 index 00000000..ba2d67dc --- /dev/null +++ b/megaradrp/processing/extractobj.py @@ -0,0 +1,256 @@ + +# Copyright 2011-2018 Universidad Complutense de Madrid +# +# This file is part of Megara DRP +# +# SPDX-License-Identifier: GPL-3.0+ +# License-Filename: LICENSE.txt +# + +"""Extract objects from RSS images""" + +import numpy +import astropy.wcs +from scipy.spatial import KDTree + + +def extract_star(rssimage, position, npoints, fiberconf, logger=None): + """Extract a star given its center and the number of fibers to extract""" + + # FIXME: handle several positions + + logger.info('extracting star') + + # fiberconf = datamodel.get_fiberconf(rssimage) + logger.debug("Configuration UUID is %s", fiberconf.conf_id) + + rssdata = rssimage[0].data + pdata = rssimage['wlmap'].data + + points = [position] + fibers = fiberconf.conected_fibers(valid_only=True) + grid_coords = [] + for fiber in fibers: + grid_coords.append((fiber.x, fiber.y)) + # setup kdtree for searching + kdtree = KDTree(grid_coords) + + # Other posibility is + # query using radius instead + # radius = 1.2 + # kdtree.query_ball_point(points, k=7, r=radius) + + dis_p, idx_p = kdtree.query(points, k=npoints) + + logger.info('Using %d nearest fibers', npoints) + totals = [] + for diss, idxs, point in zip(dis_p, idx_p, points): + # For each point + logger.info('For point %s', point) + colids = [] + coords = [] + for dis, idx in zip(diss, idxs): + fiber = fibers[idx] + colids.append(fiber.fibid - 1) + logger.debug('adding fibid %d', fiber.fibid) + + coords.append((fiber.x, fiber.y)) + + colids.sort() + flux_fiber = rssdata[colids] + flux_total = rssdata[colids].sum(axis=0) + coverage_total = pdata[colids].sum(axis=0) + + max_cover = coverage_total.max() + some_value_region = coverage_total > 0 + max_value_region = coverage_total == max_cover + valid_region = max_value_region + + # Interval with maximum coverage + nz_max, = numpy.nonzero(numpy.diff(max_value_region)) + # Interval with at least 1 fiber + nz_some, = numpy.nonzero(numpy.diff(some_value_region)) + + # Collapse the flux in the optimal region + perf = flux_fiber[:, nz_max[0] + 1: nz_max[1] + 1].sum(axis=1) + # Contribution of each fiber to the total flux, 1D + perf_norm = perf / perf.sum() + contributions = numpy.zeros(shape=(rssdata.shape[0],)) + contributions[colids] = perf_norm + # Contribution of each fiber to the total flux, 2D + flux_per_fiber = pdata * contributions[:, numpy.newaxis] + flux_sum = flux_per_fiber.sum(axis=0) + # In the region max_value_region, flux_sum == 1 + # In some_value_region 0 < flux_sum < 1 + # Outside is flux_sum == 0 + flux_correction = numpy.zeros_like(flux_sum) + flux_correction[some_value_region] = 1.0 / flux_sum[some_value_region] + + # Limit to 10 + flux_correction = numpy.clip(flux_correction, 0, 10) + + flux_total_c = flux_total * flux_correction + # plt.axvspan(nz_some[0], nz_some[1], alpha=0.2, facecolor='r') + # plt.axvspan(nz_max[0], nz_max[1], alpha=0.2, facecolor='r') + # plt.plot(flux_correction) + # plt.show() + # + # plt.axvspan(nz_some[0], nz_some[1], alpha=0.2) + # plt.axvspan(nz_max[0], nz_max[1], alpha=0.2) + # plt.plot(flux_total) + # plt.plot(flux_total * flux_correction, 'g') + # ax2 = plt.gca().twinx() + # ax2.plot(coverage_total) + # + # plt.show() + pack = flux_total_c, colids, nz_max, nz_some + # pack = flux_total_c + totals.append(pack) + + return totals[0] + + +def compute_centroid(rssdata, fiberconf, c1, c2, point, logger=None): + """Compute centroid near coordinates given by 'point'""" + + logger.debug("LCB configuration is %s", fiberconf.conf_id) + + fibers = fiberconf.conected_fibers(valid_only=True) + grid_coords = [] + for fiber in fibers: + grid_coords.append((fiber.x, fiber.y)) + # setup kdtree for searching + kdtree = KDTree(grid_coords) + + # Other posibility is + # query using radius instead + # radius = 1.2 + # kdtree.query_ball_point(points, k=7, r=radius) + + npoints = 19 + # 1 + 6 for first ring + # 1 + 6 + 12 for second ring + # 1 + 6 + 12 + 18 for third ring + points = [point] + dis_p, idx_p = kdtree.query(points, k=npoints) + + logger.info('Using %d nearest fibers', npoints) + for diss, idxs, point in zip(dis_p, idx_p, points): + # For each point + logger.info('For point %s', point) + colids = [] + coords = [] + for dis, idx in zip(diss, idxs): + fiber = fibers[idx] + colids.append(fiber.fibid - 1) + coords.append((fiber.x, fiber.y)) + + coords = numpy.asarray(coords) + flux_per_cell = rssdata[colids, c1:c2].mean(axis=1) + flux_per_cell_total = flux_per_cell.sum() + flux_per_cell_norm = flux_per_cell / flux_per_cell_total + # centroid + scf = coords.T * flux_per_cell_norm + centroid = scf.sum(axis=1) + logger.info('centroid: %s', centroid) + # central coords + c_coords = coords - centroid + scf0 = scf - centroid[:, numpy.newaxis] * flux_per_cell_norm + mc2 = numpy.dot(scf0, c_coords) + logger.info('2nd order moments, x2=%f, y2=%f, xy=%f', mc2[0,0], mc2[1,1], mc2[0,1]) + return centroid + + +def compute_dar(img, datamodel, logger=None, debug_plot=False): + """Compute Diferencial Atmospheric Refraction""" + + fiberconf = datamodel.get_fiberconf(img) + wlcalib = astropy.wcs.WCS(img[0].header) + + rssdata = img[0].data + cut1 = 500 + cut2 = 3500 + colids = [] + x = [] + y = [] + for fiber in fiberconf.fibers.values(): + colids.append(fiber.fibid - 1) + x.append(fiber.x) + y.append(fiber.y) + + cols = [] + xdar = [] + ydar = [] + delt = 50 + + point = [2.0, 2.0] + # Start in center of range + ccenter = (cut2 + cut1) // 2 + # Left + for c in range(ccenter, cut1, -delt): + c1 = c - delt // 2 + c2 = c + delt // 2 + + z = rssdata[colids, c1:c2].mean(axis=1) + centroid = compute_centroid(rssdata, fiberconf, c1, c2, point, logger=logger) + cols.append(c) + xdar.append(centroid[0]) + ydar.append(centroid[1]) + point = centroid + + cols.reverse() + xdar.reverse() + ydar.reverse() + + point = [2.0, 2.0] + # Star over + # Right + for c in range(ccenter, cut2, delt): + c1 = c - delt // 2 + c2 = c + delt // 2 + z = rssdata[colids, c1:c2].mean(axis=1) + centroid = compute_centroid(rssdata, fiberconf, c1, c2, point) + cols.append(c) + xdar.append(centroid[0]) + ydar.append(centroid[1]) + point = centroid + + rr = [[col, 0] for col in cols] + world = wlcalib.wcs_pix2world(rr, 0) + + if debug_plot: + import matplotlib.pyplot as plt + import megaradrp.visualization as vis + import megaradrp.simulation.refraction as r + from astropy import units as u + + plt.subplots_adjust(hspace=0.5) + plt.subplot(111) + ax = plt.gca() + plt.xlim([-8, 8]) + plt.ylim([-8, 8]) + col = vis.hexplot(ax, x, y, z, cmap=plt.cm.YlOrRd_r) + plt.title("Fiber map, %s %s" % (c1, c2)) + cb = plt.colorbar(col) + cb.set_label('counts') + plt.show() + + zenith_distance = 60 * u.deg + press = 79993.2 * u.Pa + rel = 0.013333333 + temp = 11.5 * u.deg_C + + ll = r.differential_p( + zenith_distance, + wl=world[:,0] * u.AA, + wl_reference=4025 * u.AA, + temperature=temp, + pressure=press, + relative_humidity=rel, + ) + plt.plot(world[:,0], xdar, '*-') + plt.plot(world[:,0], ydar, '*-') + plt.plot(world[:, 0], 2.0 * u.arcsec + ll.to(u.arcsec), '-') + plt.show() + + return world[:, 0], xdar, ydar \ No newline at end of file diff --git a/megaradrp/recipes/calibration/lcbstdstar.py b/megaradrp/recipes/calibration/lcbstdstar.py index a5c7893d..405cc9ed 100644 --- a/megaradrp/recipes/calibration/lcbstdstar.py +++ b/megaradrp/recipes/calibration/lcbstdstar.py @@ -16,12 +16,31 @@ from numina.core import Product, Parameter from numina.core.requirements import Requirement +from megaradrp.processing.extractobj import extract_star from megaradrp.recipes.scientific.base import ImageRecipe from megaradrp.types import ProcessedRSS, ProcessedFrame, ProcessedSpectrum from megaradrp.types import ReferenceSpectrumTable, ReferenceExtinctionTable from megaradrp.types import MasterSensitivity +def min_max_validator(minval=None, maxval=None): + + def checker_func(value): + import numina.exceptions + if minval is not None and value < minval: + msg = "must be >= {}".format(minval) + raise numina.exceptions.ValidationError(msg) + if maxval is not None and value > maxval: + msg = "must be <= {}".format(maxval) + raise numina.exceptions.ValidationError(msg) + return value + + return checker_func + + +nrings_check = min_max_validator(minval=1) + + class LCBStandardRecipe(ImageRecipe): """Process LCB Standard Star Recipe. @@ -61,7 +80,7 @@ class LCBStandardRecipe(ImageRecipe): """ position = Requirement(list, "Position of the reference object", default=(0, 0)) - nrings = Requirement(int, "Number of rings to extract the star", default=3) + nrings = Parameter(3, "Number of rings to extract the star", validator=nrings_check) reference_spectrum = Requirement(ReferenceSpectrumTable, "Spectrum of reference star") reference_extinction = Requirement(ReferenceExtinctionTable, "Reference extinction") sigma_resolution = Parameter(20.0, 'sigma Gaussian filter to degrade resolution ') @@ -92,10 +111,11 @@ def run(self, rinput): npoints = 1 + 3 * rinput.nrings * (rinput.nrings +1) self.logger.debug('adding %d fibers', npoints) - spectra_pack = self.extract_stars(final, rinput.position, npoints) - pack = spectra_pack[0] - # FIXME: include cover1 and cover2 - spectrum, cover1, cover2 = pack + fiberconf = self.datamodel.get_fiberconf(final) + spectra_pack = extract_star(final, rinput.position, npoints, + fiberconf, logger=self.logger) + + spectrum, colids, cover1, cover2 = spectra_pack star_spectrum = fits.PrimaryHDU(spectrum, header=final[0].header) star_interp = interp1d(rinput.reference_spectrum[:,0], rinput.reference_spectrum[:,1]) diff --git a/megaradrp/recipes/calibration/mosstdstar.py b/megaradrp/recipes/calibration/mosstdstar.py index 8b16ba8a..eed5b387 100644 --- a/megaradrp/recipes/calibration/mosstdstar.py +++ b/megaradrp/recipes/calibration/mosstdstar.py @@ -16,6 +16,7 @@ from numina.core import Product, Parameter from numina.core.requirements import Requirement +from megaradrp.processing.extractobj import extract_star from megaradrp.recipes.scientific.base import ImageRecipe from megaradrp.types import ProcessedRSS, ProcessedFrame, ProcessedSpectrum from megaradrp.types import ReferenceSpectrumTable, ReferenceExtinctionTable @@ -92,10 +93,11 @@ def run(self, rinput): npoints = 7 self.logger.debug('adding %d fibers', npoints) - spectra_pack = self.extract_stars(final, rinput.position, npoints) - pack = spectra_pack[0] - # FIXME: include cover1 and cover2 - spectrum, cover1, cover2 = pack + fiberconf = self.datamodel.get_fiberconf(final) + spectra_pack = extract_star(final, rinput.position, npoints, + fiberconf, logger=self.logger) + + spectrum, coldids, cover1, cover2 = spectra_pack star_spectrum = fits.PrimaryHDU(spectrum, header=final[0].header) star_interp = interp1d(rinput.reference_spectrum[:, 0], rinput.reference_spectrum[:, 1]) diff --git a/megaradrp/recipes/scientific/base.py b/megaradrp/recipes/scientific/base.py index d5e23f26..b16a2f7f 100644 --- a/megaradrp/recipes/scientific/base.py +++ b/megaradrp/recipes/scientific/base.py @@ -31,6 +31,7 @@ from megaradrp.processing.wavecalibration import WavelengthCalibrator from megaradrp.processing.fiberflat import Splitter, FlipLR, FiberFlatCorrector from megaradrp.processing.twilight import TwilightCorrector +from megaradrp.processing.extractobj import compute_centroid, compute_dar class ImageRecipe(MegaraBaseRecipe): @@ -113,7 +114,6 @@ def run_sky_subtraction(self, img, ignored_sky_bundles=None): sky_data[rowid] = target_data[rowid] sky_map[rowid] = target_map[rowid] if False: - plt.plot(sky_data[rowid]) plt.title("%d" % fibid) plt.show() @@ -141,341 +141,15 @@ def read_wcs(self, hdr): delt = hdr['CDELT1'] return crpix, wlr0, delt - def get_wcallib(self, lambda1, lambda2, fibras, traces, rss, neigh_info, grid): - - # Take a look at == [] - indices = [] - wlcalib = [] - for elem in traces: - if elem['aperture']['function']['coefficients']: - wlcalib.append(elem['aperture']['function']['coefficients']) - if len(indices)==0: - indices.append(0) - else: - indices.append(indices[-1]) - else: - indices.append(indices[-1]+1) - - wlcalib_aux = np.asarray(wlcalib) - final, wcsdata = self.resample_rss_flux(rss, wlcalib_aux, indices) - - hdu_f = fits.PrimaryHDU(final) - hdu_f.writeto('resample_rss.fits', clobber=True) - - fibras.sort() - - suma = np.sum(final[fibras.astype(int),lambda1:lambda2],axis=1) - sky = np.min(suma) - sumaparcial = suma - sky - - neigh_info = neigh_info[np.argsort(neigh_info[:, 0])] - - centroid_x = np.multiply(sumaparcial,neigh_info[:,2]) - centroid_x = np.sum(centroid_x, axis=0) - - centroid_y = np.multiply(sumaparcial,neigh_info[:,3]) - centroid_y = np.sum(centroid_y, axis=0) - - sumatotal = np.sum(sumaparcial, axis=0) - self.logger.info( "total sum: %s", sumatotal) - - second_order = [] - aux = np.sum(np.multiply(suma,(neigh_info[:,2] - np.mean(neigh_info[:,2]))**2),axis=0) - second_order.append(np.divide(aux ,np.sum(suma, axis=0))) - self.logger.info("Second order momentum X: %s", second_order[0]) - - aux = np.sum(np.multiply(suma,(neigh_info[:,3] - np.mean(neigh_info[:,3]))**2),axis=0) - second_order.append(np.divide(aux ,np.sum(suma, axis=0))) - self.logger.info("Second order momentum Y: %s", second_order[1]) - - aux = np.multiply(neigh_info[:,3] - np.mean(neigh_info[:,3]),neigh_info[:,2] - np.mean(neigh_info[:,2])) - aux = np.sum(np.multiply(aux,suma)) - cov = np.divide(aux ,np.sum(suma, axis=0)) - self.logger.info("Cov X,Y: %s", cov) - - centroid_x = np.divide(centroid_x, sumatotal) - self.logger.info( "centroid_x: %s", centroid_x) - - centroid_y = np.divide(centroid_y, sumatotal) - self.logger.info("centroid_y: %s", centroid_y) - - centroid = [centroid_x, centroid_y] - - peak = np.sum(final[grid.get_fiber(centroid),lambda1:lambda2],axis=0) - - return centroid, sky, peak, second_order, cov - - def generate_solution(self, points, centroid, sky, fiber, peaks, second_order, cova): - result = [] - for cont, value in enumerate(points): - lista = (value[0], value[1], centroid[cont][0],centroid[cont][1], sky[cont], fiber[cont], peaks[cont], second_order[cont][0], second_order[cont][1], cova[cont]) - result.append(lista) - return np.array(result, dtype=[('x_point','float'),('y_point','float'),('x_centroid','float'),('y_centroid','float'), ('sky','float'),('fiber','int'),('peak','float'),('x_second_order','float'), ('y_second_order','float'), ('covariance','float') ]) - - def generateJSON(self, points, centroid, sky, fiber, peaks, second_order, cova): - ''' - ''' - - self.logger.info('start JSON generation') - - result = [] - for cont, value in enumerate(points): - obj = { - 'points': value, - 'centroid': centroid[cont], - 'sky':sky[cont], - 'fiber': fiber[cont], - 'peak': peaks[cont], - 'second_order': second_order[cont], - 'covariance': cova[cont] - } - result.append(obj) - - self.logger.info('end JSON generation') - - return result - def compute_dar(self, img): import numpy.polynomial.polynomial as pol - import astropy.wcs - - fiberconf = self.datamodel.get_fiberconf(img) - wlcalib = astropy.wcs.WCS(img[0].header) - - rssdata = img[0].data - cut1 = 500 - cut2 = 3500 - colids = [] - x = [] - y = [] - for fiber in fiberconf.fibers.values(): - colids.append(fiber.fibid - 1) - x.append(fiber.x) - y.append(fiber.y) - - cols = [] - xdar = [] - ydar = [] - delt = 50 - - point = [2.0, 2.0] - # Start in center of range - ccenter = (cut2 + cut1) // 2 - # Left - for c in range(ccenter, cut1, -delt): - c1 = c - delt // 2 - c2 = c + delt // 2 - - z = rssdata[colids, c1:c2].mean(axis=1) - centroid = self.centroid(rssdata, fiberconf, c1, c2, point) - cols.append(c) - xdar.append(centroid[0]) - ydar.append(centroid[1]) - point = centroid - - cols.reverse() - xdar.reverse() - ydar.reverse() - - point = [2.0, 2.0] - # Star over - # Right - for c in range(ccenter, cut2, delt): - c1 = c - delt // 2 - c2 = c + delt // 2 - z = rssdata[colids, c1:c2].mean(axis=1) - centroid = self.centroid(rssdata, fiberconf, c1, c2, point) - cols.append(c) - xdar.append(centroid[0]) - ydar.append(centroid[1]) - point = centroid - - rr = [[col, 0] for col in cols] - world = wlcalib.wcs_pix2world(rr, 0) - - if False: - import matplotlib.pyplot as plt - import megaradrp.visualization as vis - import megaradrp.simulation.refraction as r - from astropy import units as u - - plt.subplots_adjust(hspace=0.5) - plt.subplot(111) - ax = plt.gca() - plt.xlim([-8, 8]) - plt.ylim([-8, 8]) - col = vis.hexplot(ax, x, y, z, cmap=plt.cm.YlOrRd_r) - plt.title("Fiber map, %s %s" % (c1, c2)) - cb = plt.colorbar(col) - cb.set_label('counts') - plt.show() - - zenith_distance = 60 * u.deg - press = 79993.2 * u.Pa - rel = 0.013333333 - temp = 11.5 * u.deg_C - - ll = r.differential_p( - zenith_distance, - wl=world[:,0] * u.AA, - wl_reference=4025 * u.AA, - temperature=temp, - pressure=press, - relative_humidity=rel, - ) - plt.plot(world[:,0], xdar, '*-') - plt.plot(world[:,0], ydar, '*-') - plt.plot(world[:, 0], 2.0 * u.arcsec + ll.to(u.arcsec), '-') - plt.show() - - # fit something - - print('DAR, x:', pol.polyfit(world[:, 0], xdar, deg=3)) - print('DAR: y:', pol.polyfit(world[:, 0], ydar, deg=3)) - def centroid(self, rssdata, fiberconf, c1, c2, point): - from scipy.spatial import KDTree - - self.logger.debug("LCB configuration is %s", fiberconf.conf_id) - - fibers = fiberconf.conected_fibers(valid_only=True) - grid_coords = [] - for fiber in fibers: - grid_coords.append((fiber.x, fiber.y)) - # setup kdtree for searching - kdtree = KDTree(grid_coords) - - # Other posibility is - # query using radius instead - # radius = 1.2 - # kdtree.query_ball_point(points, k=7, r=radius) - - npoints = 19 - # 1 + 6 for first ring - # 1 + 6 + 12 for second ring - # 1 + 6 + 12 + 18 for third ring - points = [point] - dis_p, idx_p = kdtree.query(points, k=npoints) - - self.logger.info('Using %d nearest fibers', npoints) - for diss, idxs, point in zip(dis_p, idx_p, points): - # For each point - self.logger.info('For point %s', point) - colids = [] - coords = [] - for dis, idx in zip(diss, idxs): - fiber = fibers[idx] - colids.append(fiber.fibid - 1) - coords.append((fiber.x, fiber.y)) - - coords = np.asarray(coords) - flux_per_cell = rssdata[colids, c1:c2].mean(axis=1) - flux_per_cell_total = flux_per_cell.sum() - flux_per_cell_norm = flux_per_cell / flux_per_cell_total - # centroid - scf = coords.T * flux_per_cell_norm - centroid = scf.sum(axis=1) - self.logger.info('centroid: %s', centroid) - # central coords - c_coords = coords - centroid - scf0 = scf - centroid[:, np.newaxis] * flux_per_cell_norm - mc2 = np.dot(scf0, c_coords) - self.logger.info('2nd order moments, x2=%f, y2=%f, xy=%f', mc2[0,0], mc2[1,1], mc2[0,1]) - return centroid - - def extract_stars(self, final, position, npoints): - from scipy.spatial import KDTree + wl, xdar, ydar = compute_dar(img, self.datamodel, logger=self.logger) + print('DAR, x:', pol.polyfit(wl, xdar, deg=3)) + print('DAR: y:', pol.polyfit(wl, ydar, deg=3)) - import matplotlib.pyplot as plt - - self.logger.info('extracting star') - - fiberconf = self.datamodel.get_fiberconf(final) - self.logger.debug("Configuration UUID is %s", fiberconf.conf_id) - rssdata = final[0].data - pdata = final['wlmap'].data - - points = [position] - fibers = fiberconf.conected_fibers(valid_only=True) - grid_coords = [] - for fiber in fibers: - grid_coords.append((fiber.x, fiber.y)) - # setup kdtree for searching - kdtree = KDTree(grid_coords) - - # Other posibility is - # query using radius instead - # radius = 1.2 - # kdtree.query_ball_point(points, k=7, r=radius) - - dis_p, idx_p = kdtree.query(points, k=npoints) - - self.logger.info('Using %d nearest fibers', npoints) - totals = [] - for diss, idxs, point in zip(dis_p, idx_p, points): - # For each point - self.logger.info('For point %s', point) - colids = [] - coords = [] - for dis, idx in zip(diss, idxs): - fiber = fibers[idx] - colids.append(fiber.fibid - 1) - self.logger.debug('adding fibid %d', fiber.fibid) - - coords.append((fiber.x, fiber.y)) - - colids.sort() - flux_fiber = rssdata[colids] - flux_total = rssdata[colids].sum(axis=0) - coverage_total = pdata[colids].sum(axis=0) - - max_cover = coverage_total.max() - some_value_region = coverage_total > 0 - max_value_region = coverage_total == max_cover - valid_region = max_value_region - - # Interval with maximum coverage - nz_max, = numpy.nonzero(numpy.diff(max_value_region)) - # Interval with at least 1 fiber - nz_some, = numpy.nonzero(numpy.diff(some_value_region)) - - # Collapse the flux in the optimal region - perf = flux_fiber[:, nz_max[0]+ 1: nz_max[1] + 1].sum(axis=1) - # Contribution of each fiber to the total flux, 1D - perf_norm = perf / perf.sum() - contributions = numpy.zeros(shape=(rssdata.shape[0],)) - contributions[colids] = perf_norm - # Contribution of each fiber to the total flux, 2D - flux_per_fiber = pdata * contributions[:, numpy.newaxis] - flux_sum = flux_per_fiber.sum(axis=0) - # In the region max_value_region, flux_sum == 1 - # In some_value_region 0 < flux_sum < 1 - # Outside is flux_sum == 0 - flux_correction = numpy.zeros_like(flux_sum) - flux_correction[some_value_region] = 1.0 / flux_sum[some_value_region] - - # Limit to 10 - flux_correction = numpy.clip(flux_correction, 0, 10) - - flux_total_c = flux_total * flux_correction - # plt.axvspan(nz_some[0], nz_some[1], alpha=0.2, facecolor='r') - # plt.axvspan(nz_max[0], nz_max[1], alpha=0.2, facecolor='r') - # plt.plot(flux_correction) - # plt.show() - # - # plt.axvspan(nz_some[0], nz_some[1], alpha=0.2) - # plt.axvspan(nz_max[0], nz_max[1], alpha=0.2) - # plt.plot(flux_total) - # plt.plot(flux_total * flux_correction, 'g') - # ax2 = plt.gca().twinx() - # ax2.plot(coverage_total) - # - # plt.show() - pack = flux_total_c, nz_max, nz_some - # pack = flux_total_c - totals.append(pack) - - return totals + def centroid(self, rssdata, fiberconf, c1, c2, point): + return compute_centroid(rssdata, fiberconf, c1, c2, point, logger=self.logger) def generate_sensitivity(self, final, spectrum, star_interp, extinc_interp, cover1, cover2, sigma=20.0):