From 7d49e6b0c573dda5e7a95c785b60f46de9698709 Mon Sep 17 00:00:00 2001 From: minzastro Date: Thu, 9 Feb 2023 23:54:03 +0100 Subject: [PATCH 1/5] Added min and median combine_functions --- reproject/mosaicking/coadd.py | 87 ++++++++++++++++++++++++++++------- 1 file changed, 70 insertions(+), 17 deletions(-) diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py index be6573c5f..b31698c15 100644 --- a/reproject/mosaicking/coadd.py +++ b/reproject/mosaicking/coadd.py @@ -6,7 +6,7 @@ from .background import determine_offset_matrix, solve_corrections_sgd from .subset_array import ReprojectedArraySubset -__all__ = ["reproject_and_coadd"] +__all__ = ['reproject_and_coadd'] def reproject_and_coadd( @@ -17,7 +17,7 @@ def reproject_and_coadd( hdu_in=None, reproject_function=None, hdu_weights=None, - combine_function="mean", + combine_function='mean', match_background=False, background_reference=None, **kwargs, @@ -92,17 +92,18 @@ def reproject_and_coadd( # Validate inputs - if combine_function not in ("mean", "sum", "median"): - raise ValueError("combine_function should be one of mean/sum/median") + if combine_function not in ('mean', 'sum', 'median', 'min'): + raise ValueError('combine_function should be one of mean/sum/median') if reproject_function is None: raise ValueError( - "reprojection function should be specified with the reproject_function argument" + 'reprojection function should be specified with the reproject_function argument' ) # Parse the output projection to avoid having to do it for each - wcs_out, shape_out = parse_output_projection(output_projection, shape_out=shape_out) + wcs_out, shape_out = parse_output_projection(output_projection, + shape_out=shape_out) # Start off by reprojecting individual images to the final projection @@ -118,7 +119,8 @@ def reproject_and_coadd( if input_weights is None: weights_in = None else: - weights_in = parse_input_weights(input_weights[idata], hdu_weights=hdu_weights) + weights_in = parse_input_weights(input_weights[idata], + hdu_weights=hdu_weights) if np.any(np.isnan(weights_in)): weights_in = np.nan_to_num(weights_in) @@ -190,7 +192,8 @@ def reproject_and_coadd( weights[reset] = 0.0 footprint *= weights - array = ReprojectedArraySubset(array, footprint, imin, imax, jmin, jmax) + array = ReprojectedArraySubset(array, footprint, + imin, imax, jmin, jmax) # TODO: make sure we gracefully handle the case where the # output image is empty (due e.g. to no overlap). @@ -213,7 +216,7 @@ def reproject_and_coadd( final_array = np.zeros(shape_out) final_footprint = np.zeros(shape_out) - if combine_function in ("mean", "sum"): + if combine_function in ('mean', 'sum'): for array in arrays: @@ -222,18 +225,68 @@ def reproject_and_coadd( # means/sums. array.array[array.footprint == 0] = 0 - final_array[array.view_in_original_array] += array.array * array.footprint + final_array[array.view_in_original_array] += \ + array.array * array.footprint final_footprint[array.view_in_original_array] += array.footprint - if combine_function == "mean": - with np.errstate(invalid="ignore"): + if combine_function == 'mean': + with np.errstate(invalid='ignore'): final_array /= final_footprint + elif combine_function == 'min': + final_array = np.ones(shape_out) * np.inf + for array in arrays: + array.array[array.footprint == 0] = np.nan + final_array[array.view_in_original_array] = \ + np.fmin(final_array[array.view_in_original_array], + array.array * array.footprint) + final_footprint[array.view_in_original_array] += array.footprint + elif combine_function == 'median': + BLOCK_SIZE = 100 + for i in np.arange(0, shape_out[0], BLOCK_SIZE): + for j in np.arange(0, shape_out[1], BLOCK_SIZE): + blocks = [] + foots = [] + for array in arrays: + part, foot = _get_block(array, i, j, BLOCK_SIZE) + part[foot == 0] = np.nan + foots.append(foot) + blocks.append(part) + final_array[i:i+BLOCK_SIZE, j:j+BLOCK_SIZE] = \ + np.nanmedian(blocks, axis=0) + final_footprint[i:i+BLOCK_SIZE, j:j+BLOCK_SIZE] = \ + np.nansum(foots, axis=0) + return final_array, final_footprint - elif combine_function == "median": - # Here we need to operate in chunks since we could otherwise run - # into memory issues +def _get_block(array, block_index_x, block_index_y, block_size=100): + """Get square block of proper dimensions from an array. - raise NotImplementedError("combine_function='median' is not yet implemented") + Args: + array (ReprojectedArraySubset): Input array to select data from. + block_index_x (int): X of the lower left corner + block_index_y (int): Y of the lower left corner + block_size (int, optional): block size in X and Y. Defaults to 100. - return final_array, final_footprint + Returns: + np.ndarray, np.ndarray: block and its integer-valued footprint. + """ + result = np.ones((block_size, block_size)) * np.nan + result_footprint = np.zeros((block_size, block_size)) + jlower = block_index_x + jupper = block_index_x + block_size + ilower = block_index_y + iupper = block_index_y + block_size + if (jlower <= array.jmax and jupper >= array.jmin and + ilower <= array.imax and iupper >= array.imin): + jlower1 = max(0, array.jmin - jlower) + jupper1 = min(jupper - jlower, array.jmax - jlower) + ilower1 = max(0, array.imin - ilower) + iupper1 = min(iupper - ilower, array.imax - ilower) + local = (slice(jlower1, jupper1), slice(ilower1, iupper1)) + external = (slice(jlower1 + jlower - array.jmin, + jupper1 + jlower - array.jmin), + slice(ilower1 + ilower - array.imin, + iupper1 + ilower - array.imin)) + result[local] = array.array[external] + result_footprint[local] = array.footprint[external] + return result, result_footprint From 23c3609fc5b1837b48fb95a413cb8ee9e1782e79 Mon Sep 17 00:00:00 2001 From: minzastro Date: Thu, 16 Feb 2023 09:18:07 +0100 Subject: [PATCH 2/5] Tests added --- reproject/mosaicking/tests/test_coadd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reproject/mosaicking/tests/test_coadd.py b/reproject/mosaicking/tests/test_coadd.py index 95e552fc6..79f36f212 100644 --- a/reproject/mosaicking/tests/test_coadd.py +++ b/reproject/mosaicking/tests/test_coadd.py @@ -78,7 +78,7 @@ def _overlapping_views(self): return views - @pytest.mark.parametrize("combine_function", ["mean", "sum"]) + @pytest.mark.parametrize("combine_function", ["mean", "sum", "min", "median"]) def test_coadd_no_overlap(self, combine_function, reproject_function): # Make sure that if all tiles are exactly non-overlapping, and From e73a57f9ee0baba92d282ec0c457bfd4ca41fafb Mon Sep 17 00:00:00 2001 From: minzastro Date: Thu, 16 Feb 2023 10:01:58 +0100 Subject: [PATCH 3/5] Tests fixed --- reproject/mosaicking/coadd.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py index b31698c15..70894e68c 100644 --- a/reproject/mosaicking/coadd.py +++ b/reproject/mosaicking/coadd.py @@ -6,6 +6,8 @@ from .background import determine_offset_matrix, solve_corrections_sgd from .subset_array import ReprojectedArraySubset +import warnings + __all__ = ['reproject_and_coadd'] @@ -251,10 +253,17 @@ def reproject_and_coadd( part[foot == 0] = np.nan foots.append(foot) blocks.append(part) - final_array[i:i+BLOCK_SIZE, j:j+BLOCK_SIZE] = \ - np.nanmedian(blocks, axis=0) - final_footprint[i:i+BLOCK_SIZE, j:j+BLOCK_SIZE] = \ - np.nansum(foots, axis=0) + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=RuntimeWarning, + message='All-NaN slice encountered') + patch_shape = final_array[i:i+BLOCK_SIZE, + j:j+BLOCK_SIZE].shape + final_array[i:i+BLOCK_SIZE, j:j+BLOCK_SIZE] = \ + np.nanmedian(blocks, axis=0)[:patch_shape[0], + :patch_shape[1]] + final_footprint[i:i+BLOCK_SIZE, j:j+BLOCK_SIZE] = \ + np.nansum(foots, axis=0)[:patch_shape[0], + :patch_shape[1]] return final_array, final_footprint From 05a9646e9f5c35a088e8abd707186dff0e8a05d3 Mon Sep 17 00:00:00 2001 From: minzastro Date: Mon, 20 Feb 2023 13:59:34 +0100 Subject: [PATCH 4/5] Comment fix --- reproject/mosaicking/coadd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py index 70894e68c..611414639 100644 --- a/reproject/mosaicking/coadd.py +++ b/reproject/mosaicking/coadd.py @@ -73,7 +73,7 @@ def reproject_and_coadd( `~astropy.io.fits.HDUList` instance, specifies the HDU to use. reproject_function : callable The function to use for the reprojection - combine_function : { 'mean', 'sum', 'median' } + combine_function : { 'mean', 'sum', 'median', 'min' } The type of function to use for combining the values into the final image. match_background : bool From 298424c4fb22916233b1e3fdaba9fae4213e34d1 Mon Sep 17 00:00:00 2001 From: minzastro Date: Wed, 3 May 2023 11:19:12 +0200 Subject: [PATCH 5/5] Fixed codestyle --- reproject/mosaicking/coadd.py | 78 +++++++++++++++++------------------ 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py index 611414639..5a08cf7a8 100644 --- a/reproject/mosaicking/coadd.py +++ b/reproject/mosaicking/coadd.py @@ -1,14 +1,14 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +import warnings + import numpy as np from ..utils import parse_input_data, parse_input_weights, parse_output_projection from .background import determine_offset_matrix, solve_corrections_sgd from .subset_array import ReprojectedArraySubset -import warnings - -__all__ = ['reproject_and_coadd'] +__all__ = ["reproject_and_coadd"] def reproject_and_coadd( @@ -19,7 +19,7 @@ def reproject_and_coadd( hdu_in=None, reproject_function=None, hdu_weights=None, - combine_function='mean', + combine_function="mean", match_background=False, background_reference=None, **kwargs, @@ -94,18 +94,17 @@ def reproject_and_coadd( # Validate inputs - if combine_function not in ('mean', 'sum', 'median', 'min'): - raise ValueError('combine_function should be one of mean/sum/median') + if combine_function not in ("mean", "sum", "median", "min"): + raise ValueError("combine_function should be one of mean/sum/median") if reproject_function is None: raise ValueError( - 'reprojection function should be specified with the reproject_function argument' + "reprojection function should be specified with the reproject_function argument" ) # Parse the output projection to avoid having to do it for each - wcs_out, shape_out = parse_output_projection(output_projection, - shape_out=shape_out) + wcs_out, shape_out = parse_output_projection(output_projection, shape_out=shape_out) # Start off by reprojecting individual images to the final projection @@ -121,8 +120,7 @@ def reproject_and_coadd( if input_weights is None: weights_in = None else: - weights_in = parse_input_weights(input_weights[idata], - hdu_weights=hdu_weights) + weights_in = parse_input_weights(input_weights[idata], hdu_weights=hdu_weights) if np.any(np.isnan(weights_in)): weights_in = np.nan_to_num(weights_in) @@ -194,8 +192,7 @@ def reproject_and_coadd( weights[reset] = 0.0 footprint *= weights - array = ReprojectedArraySubset(array, footprint, - imin, imax, jmin, jmax) + array = ReprojectedArraySubset(array, footprint, imin, imax, jmin, jmax) # TODO: make sure we gracefully handle the case where the # output image is empty (due e.g. to no overlap). @@ -218,7 +215,7 @@ def reproject_and_coadd( final_array = np.zeros(shape_out) final_footprint = np.zeros(shape_out) - if combine_function in ('mean', 'sum'): + if combine_function in ("mean", "sum"): for array in arrays: @@ -227,22 +224,21 @@ def reproject_and_coadd( # means/sums. array.array[array.footprint == 0] = 0 - final_array[array.view_in_original_array] += \ - array.array * array.footprint + final_array[array.view_in_original_array] += array.array * array.footprint final_footprint[array.view_in_original_array] += array.footprint - if combine_function == 'mean': - with np.errstate(invalid='ignore'): + if combine_function == "mean": + with np.errstate(invalid="ignore"): final_array /= final_footprint - elif combine_function == 'min': + elif combine_function == "min": final_array = np.ones(shape_out) * np.inf for array in arrays: array.array[array.footprint == 0] = np.nan - final_array[array.view_in_original_array] = \ - np.fmin(final_array[array.view_in_original_array], - array.array * array.footprint) + final_array[array.view_in_original_array] = np.fmin( + final_array[array.view_in_original_array], array.array * array.footprint + ) final_footprint[array.view_in_original_array] += array.footprint - elif combine_function == 'median': + elif combine_function == "median": BLOCK_SIZE = 100 for i in np.arange(0, shape_out[0], BLOCK_SIZE): for j in np.arange(0, shape_out[1], BLOCK_SIZE): @@ -254,16 +250,16 @@ def reproject_and_coadd( foots.append(foot) blocks.append(part) with warnings.catch_warnings(): - warnings.filterwarnings("ignore", category=RuntimeWarning, - message='All-NaN slice encountered') - patch_shape = final_array[i:i+BLOCK_SIZE, - j:j+BLOCK_SIZE].shape - final_array[i:i+BLOCK_SIZE, j:j+BLOCK_SIZE] = \ - np.nanmedian(blocks, axis=0)[:patch_shape[0], - :patch_shape[1]] - final_footprint[i:i+BLOCK_SIZE, j:j+BLOCK_SIZE] = \ - np.nansum(foots, axis=0)[:patch_shape[0], - :patch_shape[1]] + warnings.filterwarnings( + "ignore", category=RuntimeWarning, message="All-NaN slice encountered" + ) + patch_shape = final_array[i : i + BLOCK_SIZE, j : j + BLOCK_SIZE].shape + final_array[i : i + BLOCK_SIZE, j : j + BLOCK_SIZE] = np.nanmedian( + blocks, axis=0 + )[: patch_shape[0], : patch_shape[1]] + final_footprint[i : i + BLOCK_SIZE, j : j + BLOCK_SIZE] = np.nansum( + foots, axis=0 + )[: patch_shape[0], : patch_shape[1]] return final_array, final_footprint @@ -285,17 +281,21 @@ def _get_block(array, block_index_x, block_index_y, block_size=100): jupper = block_index_x + block_size ilower = block_index_y iupper = block_index_y + block_size - if (jlower <= array.jmax and jupper >= array.jmin and - ilower <= array.imax and iupper >= array.imin): + if ( + jlower <= array.jmax + and jupper >= array.jmin + and ilower <= array.imax + and iupper >= array.imin + ): jlower1 = max(0, array.jmin - jlower) jupper1 = min(jupper - jlower, array.jmax - jlower) ilower1 = max(0, array.imin - ilower) iupper1 = min(iupper - ilower, array.imax - ilower) local = (slice(jlower1, jupper1), slice(ilower1, iupper1)) - external = (slice(jlower1 + jlower - array.jmin, - jupper1 + jlower - array.jmin), - slice(ilower1 + ilower - array.imin, - iupper1 + ilower - array.imin)) + external = ( + slice(jlower1 + jlower - array.jmin, jupper1 + jlower - array.jmin), + slice(ilower1 + ilower - array.imin, iupper1 + ilower - array.imin), + ) result[local] = array.array[external] result_footprint[local] = array.footprint[external] return result, result_footprint