From cc4cb05e4b6125855911677af46e1a6246a707e3 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 29 Jan 2025 21:11:31 -0600 Subject: [PATCH 1/5] Initial antimeridian handling refactoring --- pyresample/geometry.py | 27 +++++++++++---------------- pyresample/kd_tree.py | 1 + pyresample/utils/proj4.py | 5 +++-- 3 files changed, 15 insertions(+), 18 deletions(-) diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 0925c3520..93bc2e833 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -1242,19 +1242,13 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None, projections and data cases not covering the north or south pole. The possible options are: - * "modify_extents": Set the X bounds to the edges of the data, but - add 360 to the right-most bound. This has the effect of making - the area coordinates continuous from the left side to the - right side. However, this means that some coordinates will be - outside the coordinate space of the projection. Although most - PROJ and pyresample functionality can handle this there may be - some edge cases. - * "modify_crs": Change the prime meridian of the projection - from 0 degrees longitude to 180 degrees longitude. This has - the effect of putting the data on a continuous coordinate + * "modify_extents": Deprecated. Use "modify_crs". + * "modify_crs": Modify the projection definition to wrap longitude + values so they are between 0 and 360 degrees (PROJ ``lon_wrap=180``). + This has the effect of putting the data on a continuous coordinate system. However, this means that comparing data resampled to - this resulting area and an area not over the anti-meridian - would be more difficult. + this resulting area and an area not using this longitude + wrapping may require extra care. * "global_extents": Ignore the bounds of the data and use -180/180 degrees as the west and east bounds of the data. This will generate a large output area, but with the benefit of keeping @@ -1265,6 +1259,10 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None, Shape parameters are ignored if the instance is created with the `optimize_projection` flag set to True. """ + if antimeridian_mode == "modify_extents": + warnings.warn("Antimeridian mode 'modify_extents' is deprecated. Use 'modify_crs' instead.", stacklevel=2) + antimeridian_mode = "modify_crs" + with ignore_pyproj_proj_warnings(): proj_dict = self._get_proj_dict() projection = self._projection @@ -1311,7 +1309,7 @@ def _compute_bound_centers(self, proj_dict, lonslats, antimeridian_mode): # cross anti-meridian of projection xmin, xmax = self._compute_new_x_corners_for_antimeridian(xarr, antimeridian_mode) if antimeridian_mode == "modify_crs": - proj_dict.update({"pm": 180.0}) + proj_dict.update({"lon_wrap": 180}) return proj_dict, (xmin, ymin, xmax, ymax) @staticmethod @@ -1334,9 +1332,6 @@ def _compute_new_x_corners_for_antimeridian(self, xarr, antimeridian_mode): xmax = np.nanmax(wrapped_array) if hasattr(wrapped_array, "compute"): xmin, xmax = da.compute(xmin, xmax) - if antimeridian_mode == "modify_crs": - xmin -= 180 - xmax -= 180 return xmin, xmax diff --git a/pyresample/kd_tree.py b/pyresample/kd_tree.py index 64427546e..34e55e96d 100644 --- a/pyresample/kd_tree.py +++ b/pyresample/kd_tree.py @@ -973,6 +973,7 @@ def _create_resample_kdtree(self, chunks=CHUNK_SIZE): chunks=chunks) valid_input_idx = ((source_lons >= -180) & (source_lons <= 180) & (source_lats <= 90) & (source_lats >= -90)) input_coords = lonlat2xyz(source_lons, source_lats) + input_coords = input_coords[valid_input_idx.ravel(), :] # Build kd-tree on input diff --git a/pyresample/utils/proj4.py b/pyresample/utils/proj4.py index 565643c6e..c8abfc0d5 100644 --- a/pyresample/utils/proj4.py +++ b/pyresample/utils/proj4.py @@ -180,10 +180,11 @@ def ignore_pyproj_proj_warnings(): def get_geodetic_crs_with_no_datum_shift(crs: CRS) -> CRS: """Get the geodetic CRS for the provided CRS but with no prime meridian shift.""" gcrs = crs.geodetic_crs - if gcrs.prime_meridian.longitude == 0: - return gcrs + # if gcrs.prime_meridian.longitude == 0 and gcrs.lon: + # return gcrs with ignore_pyproj_proj_warnings(): gcrs_dict = gcrs.to_dict() gcrs_dict.pop("pm", None) + # gcrs_dict.pop("lon_wrap", None) gcrs_pm0 = CRS.from_dict(gcrs_dict) return gcrs_pm0 From bd26af04350b32fc25f98fd59a3177dc116badcc Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Feb 2025 13:42:08 -0600 Subject: [PATCH 2/5] Revert deprecation of modify_extents --- pyresample/geometry.py | 22 ++++++++++++++++------ pyresample/test/test_geometry_legacy.py | 12 +++++++----- 2 files changed, 23 insertions(+), 11 deletions(-) diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 93bc2e833..be3d54594 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -1220,7 +1220,7 @@ def _get_crs_area_of_use(self, projection): return aou def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None, - antimeridian_mode=None): + antimeridian_mode: str = "modify_extents") -> AreaDefinition: """Create an AreaDefinition from this area with help of some extra info. Parameters @@ -1242,9 +1242,17 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None, projections and data cases not covering the north or south pole. The possible options are: - * "modify_extents": Deprecated. Use "modify_crs". + * "modify_extents": Set the X bounds to the edges of the data, but + add 360 to the right-most bound. This has the effect of making + the area coordinates continuous from the left side to the + right side. However, this means that some coordinates will be + outside the coordinate space of the projection. Although most + PROJ and pyresample functionality can handle this there may be + some edge cases. This is the default mode. * "modify_crs": Modify the projection definition to wrap longitude values so they are between 0 and 360 degrees (PROJ ``lon_wrap=180``). + This mode also includes the extent modifications of the + "modify_extents" mode. This has the effect of putting the data on a continuous coordinate system. However, this means that comparing data resampled to this resulting area and an area not using this longitude @@ -1258,11 +1266,13 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None, Shape parameters are ignored if the instance is created with the `optimize_projection` flag set to True. - """ - if antimeridian_mode == "modify_extents": - warnings.warn("Antimeridian mode 'modify_extents' is deprecated. Use 'modify_crs' instead.", stacklevel=2) - antimeridian_mode = "modify_crs" + .. versionchanged:: 1.33 + The "modify_crs" mode was changed to use "+lon_wrap=180" instead of + "+pm=180" as the CRS modification. This means extents of this mode + are now equivalent to the extents in "modify_extents". + + """ with ignore_pyproj_proj_warnings(): proj_dict = self._get_proj_dict() projection = self._projection diff --git a/pyresample/test/test_geometry_legacy.py b/pyresample/test/test_geometry_legacy.py index 440bd12ab..3257f675a 100644 --- a/pyresample/test/test_geometry_legacy.py +++ b/pyresample/test/test_geometry_legacy.py @@ -33,6 +33,7 @@ from pyresample import geometry from pyresample.geometry import IncompatibleAreas, combine_area_extents_vertical, concatenate_area_defs from pyresample.test.utils import catch_warnings +from pyresample.utils.proj4 import ignore_pyproj_proj_warnings class TestBaseDefinition: @@ -453,10 +454,10 @@ def test_compute_domain(self): "exclude_proj_components" ), [ - (None, (22, 60), (164.5, 24.5, 194.5, 35.5), tuple(), ("+pm=180",)), - ("modify_extents", (22, 60), (164.5, 24.5, 194.5, 35.5), tuple(), ("+pm=180",)), - ("modify_crs", (22, 60), (164.5 - 180.0, 24.5, 194.5 - 180.0, 35.5), ("+pm=180",), tuple()), - ("global_extents", (22, 720), (-180.0, 24.5, 180.0, 35.5), tuple(), ("+pm=180",)), + (None, (22, 60), (164.5, 24.5, 194.5, 35.5), tuple(), ("+lon_wrap=180", "+pm=180",)), + ("modify_extents", (22, 60), (164.5, 24.5, 194.5, 35.5), tuple(), ("+lon_wrap=180", "+pm=180")), + ("modify_crs", (22, 60), (164.5, 24.5, 194.5, 35.5), ("+lon_wrap=180",), ("+pm=180",)), + ("global_extents", (22, 720), (-180.0, 24.5, 180.0, 35.5), tuple(), ("+lon_wrap=180", "+pm=180")), ], ) @pytest.mark.parametrize("use_dask", [False, True]) @@ -471,7 +472,8 @@ def test_antimeridian_mode(self, dyn_area = geometry.DynamicAreaDefinition('test_area', '', {'proj': 'longlat'}) lons, lats = _get_fake_antimeridian_lonlats(use_dask) area = dyn_area.freeze(lonslats=(lons, lats), resolution=0.5, antimeridian_mode=antimeridian_mode) - proj_str = area.crs.to_proj4() + with ignore_pyproj_proj_warnings(): + proj_str = area.crs.to_proj4() assert area.shape == expected_shape np.testing.assert_allclose(area.area_extent, expected_extents) From b18d6c40d6c9094b2aa6715064730be458f5f97b Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 28 Feb 2025 13:44:41 -0600 Subject: [PATCH 3/5] Remove valid input/output masks in kdtree resamplers --- pyresample/future/resamplers/nearest.py | 73 ++++++++----------------- pyresample/kd_tree.py | 47 ++++++---------- pyresample/test/test_kd_tree.py | 29 +++------- 3 files changed, 49 insertions(+), 100 deletions(-) diff --git a/pyresample/future/resamplers/nearest.py b/pyresample/future/resamplers/nearest.py index e48c8b3bb..87959fe67 100644 --- a/pyresample/future/resamplers/nearest.py +++ b/pyresample/future/resamplers/nearest.py @@ -41,8 +41,8 @@ dask = None -def query_no_distance(target_lons, target_lats, valid_output_index, - mask=None, valid_input_index=None, +def query_no_distance(target_lons, target_lats, + mask=None, neighbours=None, epsilon=None, radius=None, kdtree=None): """Query the kdtree. No distances are returned. @@ -50,15 +50,11 @@ def query_no_distance(target_lons, target_lats, valid_output_index, NOTE: Dask array arguments must always come before other keyword arguments for `da.blockwise` arguments to work. """ - voi = valid_output_index - shape = voi.shape + (neighbours,) - voir = voi.ravel() if mask is not None: - mask = mask.ravel()[valid_input_index.ravel()] - target_lons_valid = target_lons.ravel()[voir] - target_lats_valid = target_lats.ravel()[voir] + mask = mask.ravel() - coords = lonlat2xyz(target_lons_valid, target_lats_valid) + # TODO: Convert between CRSes + coords = lonlat2xyz(target_lons.ravel(), target_lats.ravel()) distance_array, index_array = kdtree.query( coords, k=neighbours, @@ -72,28 +68,21 @@ def query_no_distance(target_lons, target_lats, valid_output_index, # KDTree query returns out-of-bounds neighbors as `len(arr)` # which is an invalid index, we mask those out so -1 represents # invalid values - # voi is 2D (trows, tcols) # index_array is 2D (valid output pixels, neighbors) - # there are as many Trues in voi as rows in index_array - good_pixels = index_array < kdtree.n - res_ia = np.empty(shape, dtype=int) - mask = np.zeros(shape, dtype=bool) - mask[voi, :] = good_pixels - res_ia[mask] = index_array[good_pixels] - res_ia[~mask] = -1 - return res_ia - - -def _my_index(index_arr, vii, data_arr, vii_slices=None, ia_slices=None, - fill_value=np.nan): + res_ia = index_array.astype(int, copy=False) # usually a copy since pykdtree creates uint32/64 indexes + res_ia[res_ia >= kdtree.n] = -1 + return res_ia.reshape(target_lons.shape + (neighbours,)) + + +def _my_index(index_arr, data_arr, ia_slices=None, fill_value=np.nan): """Wrap index logic for 'get_sample_from_neighbour_info' to be used inside dask map_blocks.""" - vii_slices = tuple( - x if x is not None else vii.ravel() for x in vii_slices) + # vii_slices = tuple( + # x if x is not None else vii.ravel() for x in vii_slices) mask_slices = tuple( x if x is not None else (index_arr == -1) for x in ia_slices) ia_slices = tuple( x if x is not None else index_arr for x in ia_slices) - res = data_arr[vii_slices][ia_slices] + res = data_arr[ia_slices] res[mask_slices] = fill_value return res @@ -156,21 +145,17 @@ def _create_resample_kdtree(self, chunks=CHUNK_SIZE): """Set up kd tree on input.""" source_lons, source_lats = self.source_geo_def.get_lonlats( chunks=chunks) - valid_input_idx = ((source_lons >= -180) & (source_lons <= 180) & (source_lats <= 90) & (source_lats >= -90)) input_coords = lonlat2xyz(source_lons, source_lats) - input_coords = input_coords[valid_input_idx.ravel(), :] # Build kd-tree on input - input_coords = input_coords.astype(np.float64) + input_coords = input_coords.astype(np.float64, copy=False) delayed_kdtree = dask.delayed(KDTree, pure=True)(input_coords) - return valid_input_idx, delayed_kdtree + return delayed_kdtree def _query_resample_kdtree(self, resample_kdtree, tlons, tlats, - valid_input_index, - valid_output_index, mask, neighbors, radius_of_influence, @@ -181,12 +166,12 @@ def _query_resample_kdtree(self, else: ndims = self.source_geo_def.ndim dims = 'mn'[:ndims] - args = (mask, dims, valid_input_index, dims) + args = (mask, dims, dims) # res.shape = rows, cols, neighbors # j=rows, i=cols, k=neighbors, m=source rows, n=source cols res = da.blockwise( query_no_distance, 'jik', tlons, 'ji', tlats, 'ji', - valid_output_index, 'ji', *args, kdtree=resample_kdtree, + *args, kdtree=resample_kdtree, neighbours=neighbors, epsilon=epsilon, radius=radius_of_influence, dtype=np.int64, meta=np.array((), dtype=np.int64), @@ -208,27 +193,25 @@ def _get_neighbor_info(self, mask, neighbors, radius_of_influence, epsilon): # Create kd-tree chunks = mask.chunks if mask is not None else CHUNK_SIZE - valid_input_idx, resample_kdtree = self._create_resample_kdtree(chunks=chunks) + resample_kdtree = self._create_resample_kdtree(chunks=chunks) # TODO: Add 'chunks' keyword argument to this method and use it target_lons, target_lats = self.target_geo_def.get_lonlats(chunks=CHUNK_SIZE) - valid_output_idx = ((target_lons >= -180) & (target_lons <= 180) & (target_lats <= 90) & (target_lats >= -90)) if mask is not None: if mask.shape != self.source_geo_def.shape: raise ValueError("'mask' must be the same shape as the source geo definition") mask = mask.data index_arr = self._query_resample_kdtree( - resample_kdtree, target_lons, target_lats, valid_input_idx, - valid_output_idx, mask, + resample_kdtree, target_lons, target_lats, + mask, neighbors, radius_of_influence, epsilon) - return valid_input_idx, index_arr + return index_arr def get_sample_from_neighbor_info( self, data, - valid_input_index, index_array, neighbors=1, fill_value=np.nan): @@ -278,7 +261,6 @@ def get_sample_from_neighbor_info( "handle more than 1 neighbor yet.") # Convert from multiple neighbor shape to 1 neighbor ia = index_array[:, :, 0] - vii = valid_input_index src_geo_dims = self._get_src_geo_dims() dst_geo_dims = ('y', 'x') @@ -287,7 +269,6 @@ def get_sample_from_neighbor_info( # shape of the source data after we flatten the geo dimensions flat_src_shape = [] # slice objects to index in to the source data - vii_slices = [] ia_slices = [] # whether we have seen the geo dims in our analysis geo_handled = False @@ -302,7 +283,6 @@ def get_sample_from_neighbor_info( src_dim_to_ind[dim] = i if dim in src_geo_dims and not geo_handled: flat_src_shape.append(-1) - vii_slices.append(None) # mark for replacement ia_slices.append(None) # mark for replacement flat_adim.append(i) src_adims.append(i) @@ -310,7 +290,6 @@ def get_sample_from_neighbor_info( geo_handled = True elif dim not in src_geo_dims: flat_src_shape.append(data.sizes[dim]) - vii_slices.append(slice(None)) ia_slices.append(slice(None)) src_adims.append(i) dst_dims.append(dim) @@ -322,7 +301,6 @@ def get_sample_from_neighbor_info( # neighbors_dim = i + 3 new_data = data.data.reshape(flat_src_shape) - vii = vii.ravel() dst_adims = [dst_dim_to_ind[dim] for dim in dst_dims] ia_adims = [dst_dim_to_ind[dim] for dim in dst_geo_dims] # FUTURE: when we allow more than one neighbor add neighbors dimension @@ -336,9 +314,8 @@ def get_sample_from_neighbor_info( res = da.blockwise( _my_index, dst_adims, ia, ia_adims, - vii, flat_adim, new_data, src_adims, - vii_slices=vii_slices, ia_slices=ia_slices, + ia_slices=ia_slices, fill_value=fill_value, meta=np.array((), dtype=new_data.dtype), dtype=new_data.dtype, concatenate=True) @@ -410,10 +387,9 @@ def precompute(self, mask=None, radius_of_influence=None, epsilon=0): internal_cache_key = (mask_hash, neighbors, radius_of_influence, epsilon) in_int_cache = internal_cache_key in self._internal_cache if not in_int_cache: - valid_input_index, index_arr = self._get_neighbor_info( + index_arr = self._get_neighbor_info( mask, neighbors, radius_of_influence, epsilon) item_to_cache = { - "valid_input_index": valid_input_index, "index_array": index_arr, } self._internal_cache[internal_cache_key] = item_to_cache @@ -465,7 +441,6 @@ def resample(self, data, mask_area=None, fill_value=np.nan, precompute_dict = self._internal_cache[cache_key] result = self.get_sample_from_neighbor_info( new_data, - precompute_dict["valid_input_index"], precompute_dict["index_array"], fill_value=fill_value) return self._verify_result_object_type(result, data) diff --git a/pyresample/kd_tree.py b/pyresample/kd_tree.py index 34e55e96d..ea06d1500 100644 --- a/pyresample/kd_tree.py +++ b/pyresample/kd_tree.py @@ -899,7 +899,7 @@ def _prepare_result(result, output_shape, is_masked_data, use_masked_fill_value, return result -class XArrayResamplerNN(object): +class XArrayResamplerNN: """Resampler for Xarray DataArray objects with the nearest neighbor algorithm.""" def __init__(self, @@ -969,23 +969,19 @@ def _compute_radius_of_influence(self): def _create_resample_kdtree(self, chunks=CHUNK_SIZE): """Set up kd tree on input.""" - source_lons, source_lats = self.source_geo_def.get_lonlats( - chunks=chunks) - valid_input_idx = ((source_lons >= -180) & (source_lons <= 180) & (source_lats <= 90) & (source_lats >= -90)) + source_lons, source_lats = self.source_geo_def.get_lonlats(chunks=chunks) + # TODO: Pass XY coordinates of input geometry along with CRS input_coords = lonlat2xyz(source_lons, source_lats) - input_coords = input_coords[valid_input_idx.ravel(), :] - # Build kd-tree on input - input_coords = input_coords.astype(np.float64) + input_coords = input_coords.astype(np.float64, copy=False) delayed_kdtree = dask.delayed(KDTree, pure=True)(input_coords) - return valid_input_idx, delayed_kdtree + return delayed_kdtree def query_resample_kdtree(self, resample_kdtree, tlons, tlats, - valid_oi, mask): """Query kd-tree on slice of target coordinates.""" if mask is None: @@ -993,11 +989,11 @@ def query_resample_kdtree(self, else: ndims = self.source_geo_def.ndim dims = 'mn'[:ndims] - args = (mask, dims, self.valid_input_index, dims) + args = (mask, dims) # res.shape = rows, cols, neighbors # j=rows, i=cols, k=neighbors, m=source rows, n=source cols res = blockwise(query_no_distance, 'jik', tlons, 'ji', tlats, 'ji', - valid_oi, 'ji', *args, kdtree=resample_kdtree, + *args, kdtree=resample_kdtree, neighbours=self.neighbours, epsilon=self.epsilon, radius=self.radius_of_influence, dtype=np.int64, meta=np.array((), dtype=np.int64), @@ -1019,29 +1015,23 @@ def get_neighbour_info(self, mask=None): # Create kd-tree chunks = mask.chunks if mask is not None else CHUNK_SIZE - valid_input_idx, resample_kdtree = self._create_resample_kdtree( - chunks=chunks) - self.valid_input_index = valid_input_idx + resample_kdtree = self._create_resample_kdtree(chunks=chunks) self.delayed_kdtree = resample_kdtree # TODO: Add 'chunks' keyword argument to this method and use it + # TODO: Use XY coordinates and pass CRS target_lons, target_lats = self.target_geo_def.get_lonlats(chunks=CHUNK_SIZE) - valid_output_idx = ((target_lons >= -180) & (target_lons <= 180) & (target_lats <= 90) & (target_lats >= -90)) if mask is not None: if mask.shape != self.source_geo_def.shape: raise ValueError("'mask' must be the same shape as the source geo definition") mask = mask.data - index_arr, distance_arr = self.query_resample_kdtree( - resample_kdtree, target_lons, target_lats, valid_output_idx, mask) + index_arr, distance_arr = self.query_resample_kdtree(resample_kdtree, target_lons, target_lats, mask) - self.valid_output_index, self.index_array = valid_output_idx, index_arr + self.index_array = index_arr self.distance_array = distance_arr - return (self.valid_input_index, - self.valid_output_index, - self.index_array, - self.distance_array) + return self.index_array, self.distance_array def get_sample_from_neighbour_info(self, data, fill_value=np.nan): """Get the pixels matching the target area. @@ -1085,7 +1075,7 @@ def get_sample_from_neighbour_info(self, data, fill_value=np.nan): "handle more than 1 neighbor yet.") # Convert from multiple neighbor shape to 1 neighbor ia = self.index_array[:, :, 0] - vii = self.valid_input_index + # vii = self.valid_input_index src_geo_dims, dst_geo_dims = self._get_valid_dims(data) # FIXME: Can't include coordinates whose dimensions depend on the geo dims either @@ -1105,7 +1095,7 @@ def contain_coords(var, coord_list): # shape of the source data after we flatten the geo dimensions flat_src_shape = [] # slice objects to index in to the source data - vii_slices = [] + # vii_slices = [] ia_slices = [] # whether we have seen the geo dims in our analysis geo_handled = False @@ -1120,7 +1110,7 @@ def contain_coords(var, coord_list): src_dim_to_ind[dim] = i if dim in src_geo_dims and not geo_handled: flat_src_shape.append(-1) - vii_slices.append(None) # mark for replacement + # vii_slices.append(None) # mark for replacement ia_slices.append(None) # mark for replacement flat_adim.append(i) src_adims.append(i) @@ -1128,7 +1118,7 @@ def contain_coords(var, coord_list): geo_handled = True elif dim not in src_geo_dims: flat_src_shape.append(data.sizes[dim]) - vii_slices.append(slice(None)) + # vii_slices.append(slice(None)) ia_slices.append(slice(None)) src_adims.append(i) dst_dims.append(dim) @@ -1140,7 +1130,7 @@ def contain_coords(var, coord_list): # neighbors_dim = i + 3 new_data = data.data.reshape(flat_src_shape) - vii = vii.ravel() + # vii = vii.ravel() dst_adims = [dst_dim_to_ind[dim] for dim in dst_dims] ia_adims = [dst_dim_to_ind[dim] for dim in dst_geo_dims] # FUTURE: when we allow more than one neighbor add neighbors dimension @@ -1153,9 +1143,8 @@ def contain_coords(var, coord_list): # then we can avoid all of this complicated blockwise stuff res = blockwise(_my_index, dst_adims, ia, ia_adims, - vii, flat_adim, new_data, src_adims, - vii_slices=vii_slices, ia_slices=ia_slices, + ia_slices=ia_slices, fill_value=fill_value, meta=np.array((), dtype=new_data.dtype), dtype=new_data.dtype, concatenate=True) diff --git a/pyresample/test/test_kd_tree.py b/pyresample/test/test_kd_tree.py index b1e52e66a..76e35ec9d 100644 --- a/pyresample/test/test_kd_tree.py +++ b/pyresample/test/test_kd_tree.py @@ -827,9 +827,7 @@ def test_nearest_swath_1d_mask_to_grid_1n(self): neighbours=1) data = self.tdata_1d ninfo = resampler.get_neighbour_info(mask=data.isnull()) - for val in ninfo[:3]: - # vii, ia, voi - self.assertIsInstance(val, da.Array) + self.assertIsInstance(ninfo[0], da.Array) # index array res = resampler.get_sample_from_neighbour_info(data) self.assertIsInstance(res, xr.DataArray) self.assertIsInstance(res.data, da.Array) @@ -851,14 +849,11 @@ def test_nearest_type_preserve(self): resampler = XArrayResamplerNN(self.tswath_1d, self.tgrid, radius_of_influence=100000, neighbours=1) - data = self.tdata_1d data = xr.DataArray(da.from_array(np.array([1, 2, 3]), chunks=5), dims=('my_dim1',)) ninfo = resampler.get_neighbour_info() - for val in ninfo[:3]: - # vii, ia, voi - self.assertIsInstance(val, da.Array) + self.assertIsInstance(ninfo[0], da.Array) # index array res = resampler.get_sample_from_neighbour_info(data, fill_value=255) self.assertIsInstance(res, xr.DataArray) self.assertIsInstance(res.data, da.Array) @@ -883,9 +878,7 @@ def test_nearest_swath_2d_mask_to_area_1n(self): radius_of_influence=50000, neighbours=1) ninfo = resampler.get_neighbour_info(mask=data.isnull()) - for val in ninfo[:3]: - # vii, ia, voi - self.assertIsInstance(val, da.Array) + self.assertIsInstance(ninfo[0], da.Array) # index array res = resampler.get_sample_from_neighbour_info(data) self.assertIsInstance(res, xr.DataArray) self.assertIsInstance(res.data, da.Array) @@ -907,9 +900,7 @@ def test_nearest_area_2d_to_area_1n(self): neighbours=1) with assert_maximum_dask_computes(0): ninfo = resampler.get_neighbour_info() - for val in ninfo[:3]: - # vii, ia, voi - self.assertIsInstance(val, da.Array) + self.assertIsInstance(ninfo[0], da.Array) # index array with pytest.raises(ValueError): resampler.get_sample_from_neighbour_info(data) @@ -934,9 +925,7 @@ def test_nearest_area_2d_to_area_1n_no_roi(self): resampler = XArrayResamplerNN(self.src_area_2d, self.area_def, neighbours=1) ninfo = resampler.get_neighbour_info() - for val in ninfo[:3]: - # vii, ia, voi - self.assertIsInstance(val, da.Array) + self.assertIsInstance(ninfo[0], da.Array) # index array with pytest.raises(ValueError): resampler.get_sample_from_neighbour_info(data) @@ -977,9 +966,7 @@ def test_nearest_area_2d_to_area_1n_3d_data(self): radius_of_influence=50000, neighbours=1) ninfo = resampler.get_neighbour_info() - for val in ninfo[:3]: - # vii, ia, voi - self.assertIsInstance(val, da.Array) + self.assertIsInstance(ninfo[0], da.Array) # index array with pytest.raises(ValueError): resampler.get_sample_from_neighbour_info(data) @@ -1006,9 +993,7 @@ def test_nearest_swath_1d_mask_to_grid_8n(self): neighbours=8) data = self.tdata_1d ninfo = resampler.get_neighbour_info(mask=data.isnull()) - for val in ninfo[:3]: - # vii, ia, voi - self.assertIsInstance(val, da.Array) + self.assertIsInstance(ninfo[0], da.Array) # index array res = resampler.get_sample_from_neighbour_info(data) self.assertIsInstance(res, xr.DataArray) self.assertIsInstance(res.data, da.Array) From 06057aa71da6dbf5a18ccd586dfd4cedd648ec1e Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 3 Mar 2025 20:30:35 -0600 Subject: [PATCH 4/5] Add test for lon_wrap=180 in nearest resampler --- pyresample/test/conftest.py | 14 ++++++++++++++ pyresample/test/test_resamplers/test_nearest.py | 14 +++++++++----- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/pyresample/test/conftest.py b/pyresample/test/conftest.py index 0bce5f63f..24253cb15 100644 --- a/pyresample/test/conftest.py +++ b/pyresample/test/conftest.py @@ -254,6 +254,20 @@ def area_def_lonlat_pm180_target(): ) +@pytest.fixture(scope="session") +def area_def_lonlat_lonwrap180_target(): + """Create an AreaDefinition with a geographic lon/lat projection wrapping around the antimeridian (800, 850).""" + return AreaDefinition( + { + 'proj': 'longlat', + 'lon_wrap': '180.0', + 'datum': 'WGS84', + 'no_defs': None, + }, + (DST_AREA_SHAPE[0], DST_AREA_SHAPE[1]), + (160.0, 20.0, 200.0, 35.0) + ) + @pytest.fixture(scope="session") def coord_def_2d_float32_dask(): """Create a 2D CoordinateDefinition of dask arrays (4, 3).""" diff --git a/pyresample/test/test_resamplers/test_nearest.py b/pyresample/test/test_resamplers/test_nearest.py index a4a16db51..e32a08588 100644 --- a/pyresample/test/test_resamplers/test_nearest.py +++ b/pyresample/test/test_resamplers/test_nearest.py @@ -124,15 +124,19 @@ def test_nearest_area_2d_to_area_1n(self, area_def_stere_source, data_2d_float32 assert cross_sum == expected assert res.shape == resampler.target_geo_def.shape - def test_nearest_swath_2d_to_area_1n_pm180(self, swath_def_2d_xarray_dask_antimeridian, data_2d_float32_xarray_dask, - area_def_lonlat_pm180_target): + @pytest.mark.parametrize("dst_area", [lf("area_def_lonlat_pm180_target"), lf("area_def_lonlat_lonwrap180_target")]) + def test_nearest_swath_2d_to_area_1n_pm180( + self, + swath_def_2d_xarray_dask_antimeridian, + data_2d_float32_xarray_dask, + dst_area + ): """Test 2D swath definition to 2D area definition; 1 neighbor; output prime meridian at 180 degrees.""" - resampler = KDTreeNearestXarrayResampler( - swath_def_2d_xarray_dask_antimeridian, area_def_lonlat_pm180_target) + resampler = KDTreeNearestXarrayResampler(swath_def_2d_xarray_dask_antimeridian, dst_area) res = resampler.resample(data_2d_float32_xarray_dask, radius_of_influence=50000) assert isinstance(res, xr.DataArray) assert isinstance(res.data, da.Array) - _check_common_metadata(res, isinstance(area_def_lonlat_pm180_target, AreaDefinition)) + _check_common_metadata(res, isinstance(dst_area, AreaDefinition)) res = res.values cross_sum = float(np.nansum(res)) expected = 115591.0 From 2c4004bc9b75a3a58bb2280220fc0615dd6dd0a4 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 4 Mar 2025 15:49:26 -0600 Subject: [PATCH 5/5] Add antimeridian kdtree tests with antimeridian input --- pyresample/test/conftest.py | 20 +++++++++++++++++ .../test/test_resamplers/test_nearest.py | 22 ++++++++++++++++++- 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/pyresample/test/conftest.py b/pyresample/test/conftest.py index 24253cb15..198161e56 100644 --- a/pyresample/test/conftest.py +++ b/pyresample/test/conftest.py @@ -201,6 +201,26 @@ def area_def_lcc_conus_1km(): return area_def +@pytest.fixture(scope="session") +def area_def_lonlat_pm180(): + """Create an AreaDefinition with an lon/lat projection over the antimeridian with pm shift (1500, 2000).""" + proj_str = "+proj=longlat +pm=180 +datum=WGS84 +no_defs" + crs = CRS.from_string(proj_str) + area_def = AreaDefinition(crs, (SRC_AREA_SHAPE[0], SRC_AREA_SHAPE[1]), + (-10, 15, 20, 30)) + return area_def + + +@pytest.fixture(scope="session") +def area_def_lonlat_lonwrap180(): + """Create an AreaDefinition with an lon/lat projection over the antimeridian (1500, 2000).""" + proj_str = "+proj=longlat +lon_wrap=180 +datum=WGS84 +no_defs" + crs = CRS.from_string(proj_str) + area_def = AreaDefinition(crs, (SRC_AREA_SHAPE[0], SRC_AREA_SHAPE[1]), + (170, 15, 200, 30)) + return area_def + + @pytest.fixture(scope="session") def area_def_stere_source(): """Create an AreaDefinition with a polar-stereographic projection (10, 50). diff --git a/pyresample/test/test_resamplers/test_nearest.py b/pyresample/test/test_resamplers/test_nearest.py index e32a08588..af6d2062c 100644 --- a/pyresample/test/test_resamplers/test_nearest.py +++ b/pyresample/test/test_resamplers/test_nearest.py @@ -129,7 +129,7 @@ def test_nearest_swath_2d_to_area_1n_pm180( self, swath_def_2d_xarray_dask_antimeridian, data_2d_float32_xarray_dask, - dst_area + dst_area, ): """Test 2D swath definition to 2D area definition; 1 neighbor; output prime meridian at 180 degrees.""" resampler = KDTreeNearestXarrayResampler(swath_def_2d_xarray_dask_antimeridian, dst_area) @@ -143,6 +143,26 @@ def test_nearest_swath_2d_to_area_1n_pm180( assert cross_sum == expected assert res.shape == resampler.target_geo_def.shape + @pytest.mark.parametrize("src_area", [lf("area_def_lonlat_pm180"), lf("area_def_lonlat_lonwrap180")]) + @pytest.mark.parametrize("dst_area", [lf("area_def_lonlat_pm180_target"), lf("area_def_lonlat_lonwrap180_target")]) + def test_nearest_area_2d_to_area_1n_pm180( + self, + src_area, + data_2d_float32_xarray_dask, + dst_area, + ): + """Test 2D swath definition to 2D area definition; 1 neighbor; over the antimeridian.""" + resampler = KDTreeNearestXarrayResampler(src_area, dst_area) + res = resampler.resample(data_2d_float32_xarray_dask, radius_of_influence=50000) + assert isinstance(res, xr.DataArray) + assert isinstance(res.data, da.Array) + _check_common_metadata(res, isinstance(dst_area, AreaDefinition)) + res = res.values + cross_sum = float(np.nansum(res)) + expected = 78499.0 + assert cross_sum == expected + assert res.shape == resampler.target_geo_def.shape + def test_nearest_area_2d_to_area_1n_no_roi(self, area_def_stere_source, data_2d_float32_xarray_dask, area_def_stere_target): """Test 2D area definition to 2D area definition; 1 neighbor, no radius of influence."""