From 599ce9e6c4dca52db60b07587c68a9131cef45ab Mon Sep 17 00:00:00 2001 From: Samuel Scherrer Date: Thu, 26 Nov 2020 12:09:03 +0100 Subject: [PATCH 1/6] transform lon to be in (-180,180] in BasicGrid --- src/pygeogrids/grids.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/pygeogrids/grids.py b/src/pygeogrids/grids.py index 45da94f..9ff71d2 100644 --- a/src/pygeogrids/grids.py +++ b/src/pygeogrids/grids.py @@ -31,6 +31,7 @@ import numpy as np import numpy.testing as nptest +import warnings try: from osgeo import ogr ogr_installed = True @@ -93,6 +94,11 @@ class BasicGrid(object): The shape has to be given as (lat2d, lon2d) It it is not given the shape is set to the length of the input lon and lat arrays. + transform_lon : bool or None, optional (default: None) + Whether to transform longitudes to values between -180 and 180. + By default values are transformed, but a warning is issued. + To turn off the warning, set this to ``True``, to turn of + transformation set this to ``False``. Attributes ---------- @@ -148,7 +154,7 @@ class BasicGrid(object): """ def __init__(self, lon, lat, gpis=None, geodatum='WGS84', subset=None, - setup_kdTree=True, shape=None): + setup_kdTree=True, shape=None, transform_lon=None): """ init method, prepares lon and lat arrays for _transform_lonlats if necessary @@ -167,6 +173,18 @@ def __init__(self, lon, lat, gpis=None, geodatum='WGS84', subset=None, self.n_gpi = len(lon) + # transfrom longitudes to be between -180 and 180 if they are between 0 + # and 360 + if transform_lon or transform_lon is None: + if np.any(lon > 180): + lon[lon > 180] -= 360 + if transform_lon is None: + warnings.warn( + "Longitude values have been transformed to be in" + " (-180, 180]. If this was not intended or to suppress" + " this warning set the transform_lon keyword argument" + ) + self.arrlon = lon self.arrlat = lat From 0293d3e3fd4274743e12fcc9a5cbe46d1457420e Mon Sep 17 00:00:00 2001 From: Samuel Scherrer Date: Fri, 27 Nov 2020 09:00:46 +0100 Subject: [PATCH 2/6] fix #65 --- src/pygeogrids/netcdf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pygeogrids/netcdf.py b/src/pygeogrids/netcdf.py index aba3be0..6972371 100644 --- a/src/pygeogrids/netcdf.py +++ b/src/pygeogrids/netcdf.py @@ -87,8 +87,8 @@ def save_lonlat(filename, arrlon, arrlat, geodatum, arrcell=None, type(global_attrs['shape']) is not int and len(global_attrs['shape']) == 2): - latsize = global_attrs['shape'][1] - lonsize = global_attrs['shape'][0] + latsize = global_attrs['shape'][0] + lonsize = global_attrs['shape'][1] ncfile.createDimension("lat", latsize) ncfile.createDimension("lon", lonsize) gpisize = global_attrs['shape'][0] * global_attrs['shape'][1] From 7e43bbb1d502d98553f0245abda3e7e7d5bffbb6 Mon Sep 17 00:00:00 2001 From: Samuel Scherrer Date: Fri, 27 Nov 2020 09:26:33 +0100 Subject: [PATCH 3/6] added test for transform_lon --- tests/test_grid.py | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/tests/test_grid.py b/tests/test_grid.py index d50ff0e..08ad010 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -33,8 +33,10 @@ import numpy.testing as nptest import numpy as np from osgeo import ogr +import pytest +import warnings -from pygeogrids.grids import lonlat2cell +from pygeogrids.grids import lonlat2cell, BasicGrid import pygeogrids as grids @@ -659,5 +661,30 @@ def test_shpgrid(self): assert subgrid.activearrlat == 46 +@pytest.mark.filterwarnings("error") +def test_BasicGrid_transform_lon(): + """ + Tests whether transforming longitudes works as expected. + """ + + lat = np.asarray([10, -10, 5, 42]) + lon_pos = np.asarray([0, 90, 180, 270]) + lon_centered = np.asarray([0, 90, 180, -90]) + + # case 1: warning and transformation + with pytest.warns(UserWarning): + grid = BasicGrid(lon_pos, lat) + assert np.all(grid.arrlon == lon_centered) + + # case 2: no warning and transform + grid = BasicGrid(lon_pos, lat, transform_lon=True) + assert np.all(grid.arrlon == lon_centered) + + # case 3: no warning and no transform + grid = BasicGrid(lon_pos, lat, transform_lon=False) + assert np.all(grid.arrlon == lon_pos) + + + if __name__ == "__main__": unittest.main() From df8cea54bcf80bce7d6f94bf06bfa1b832abbf78 Mon Sep 17 00:00:00 2001 From: Samuel Scherrer Date: Fri, 27 Nov 2020 09:28:18 +0100 Subject: [PATCH 4/6] fix #64 --- src/pygeogrids/grids.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pygeogrids/grids.py b/src/pygeogrids/grids.py index 9ff71d2..24eabb7 100644 --- a/src/pygeogrids/grids.py +++ b/src/pygeogrids/grids.py @@ -386,7 +386,7 @@ def find_nearest_gpi(self, lon, lat, max_dist=np.Inf): lat : float or iterable Latitude of point. max_dist : float, optional - Maximum distance to consider for search (default: np.Inf). + Maximum distance [m] to consider for search (default: np.Inf). Returns ------- @@ -397,7 +397,7 @@ def find_nearest_gpi(self, lon, lat, max_dist=np.Inf): At the moment not on a great circle but in spherical cartesian coordinates. """ - gpi, distance = self.find_k_nearest_gpi(lon, lat, max_dist=np.Inf, k=1) + gpi, distance = self.find_k_nearest_gpi(lon, lat, max_dist=max_dist, k=1) if not _element_iterable(lon): gpi = gpi[0] From 26256fd8970091894467911a09855f1927646ee6 Mon Sep 17 00:00:00 2001 From: Samuel Scherrer Date: Fri, 27 Nov 2020 09:40:22 +0100 Subject: [PATCH 5/6] fixed tests according to #65 --- tests/test_netcdf.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_netcdf.py b/tests/test_netcdf.py index 7183a31..cb1bc5e 100644 --- a/tests/test_netcdf.py +++ b/tests/test_netcdf.py @@ -47,12 +47,12 @@ def setUp(self): self.subset = np.sort(np.random.choice(np.arange(self.lats.size), size=500, replace=False)) self.basic = grids.BasicGrid(self.lons, self.lats, subset=self.subset, - shape=(360, 180)) + shape=(180, 360)) self.basic_shape_gpis = grids.BasicGrid(self.lons, self.lats, gpis=np.arange(self.lats.size), subset=self.subset, - shape=(360, 180)) + shape=(180, 360)) self.basic_generated = grids.genreg_grid(1, 1) self.basic_irregular = grids.BasicGrid(np.random.random(360 * 180) * 360 - 180, np.random.random( @@ -63,7 +63,7 @@ def setUp(self): self.cellgrid_shape = grids.CellGrid(self.lons, self.lats, self.cells, subset=self.subset, - shape=(360, 180)) + shape=(180, 360)) self.testfile = tempfile.NamedTemporaryFile().name @@ -106,8 +106,8 @@ def test_save_basicgrid_generated(self): nptest.assert_array_equal(sorted(self.basic.gpis[self.subset]), sorted(nc_data.variables['gpi'][:].flatten()[stored_subset])) assert nc_data.test == 'test_attribute' - assert nc_data.shape[1] == 180 - assert nc_data.shape[0] == 360 + assert nc_data.shape[0] == 180 + assert nc_data.shape[1] == 360 def test_save_basicgrid_irregular_nc(self): grid_nc.save_grid(self.testfile, From 4679d8fbe6844ce3ec8e6e94af2af4753922de23 Mon Sep 17 00:00:00 2001 From: Samuel Scherrer Date: Fri, 27 Nov 2020 10:05:22 +0100 Subject: [PATCH 6/6] tests for max_dist in find_nearest_gpi The tests revealed that d=inf and gpi=len(gpis)+1 is returned in case no gpi is found within max_dist. I added a check for infinite distances in nearest_neighbor.find_nearest_index. --- src/pygeogrids/grids.py | 2 +- src/pygeogrids/nearest_neighbor.py | 24 +++++++++++++++++++----- tests/test_grid.py | 13 +++++++++++++ 3 files changed, 33 insertions(+), 6 deletions(-) diff --git a/src/pygeogrids/grids.py b/src/pygeogrids/grids.py index 24eabb7..3a04fa6 100644 --- a/src/pygeogrids/grids.py +++ b/src/pygeogrids/grids.py @@ -399,7 +399,7 @@ def find_nearest_gpi(self, lon, lat, max_dist=np.Inf): """ gpi, distance = self.find_k_nearest_gpi(lon, lat, max_dist=max_dist, k=1) - if not _element_iterable(lon): + if not _element_iterable(lon) and len(gpi) > 0: gpi = gpi[0] distance = distance[0] diff --git a/src/pygeogrids/nearest_neighbor.py b/src/pygeogrids/nearest_neighbor.py index e088e95..45e4c3d 100644 --- a/src/pygeogrids/nearest_neighbor.py +++ b/src/pygeogrids/nearest_neighbor.py @@ -178,12 +178,22 @@ def find_nearest_index(self, lon, lat, max_dist=np.Inf, k=1): great circle distance at the moment. This should be OK for most applications that look for the nearest neighbor which should not be hundreds of kilometers away. + If no point was found within the maximum distance to consider, an + empty array is returned. ind : int, numpy.array - indices of nearest neighbor + If ``self.grid`` is ``False`` indices of nearest neighbor. + If no point was found within the maximum distance to consider, an + empty array is returned. index_lon : numpy.array, optional - if self.grid is True then return index into lon array of grid definition + If ``self.grid`` is ``True`` then return index into lon array of + grid definition. + If no point was found within the maximum distance to consider, an + empty array is returned. index_lat : numpy.array, optional - if self.grid is True then return index into lat array of grid definition + If ``self.grid`` is ``True`` then return index into lat array of + grid definition. + If no point was found within the maximum distance to consider, an + empty array is returned. """ if self.kdtree is None: self._build_kdtree() @@ -193,11 +203,15 @@ def find_nearest_index(self, lon, lat, max_dist=np.Inf, k=1): d, ind = self.kdtree.query( query_coords, distance_upper_bound=max_dist, k=k) + # if no point was found, d == inf + if not np.all(np.isfinite(d)): + d, ind = np.array([]), np.array([]) + if not self.grid: return d, ind else: - # calculate index position in grid definition arrays assuming row-major - # flattening of arrays after numpy.meshgrid + # calculate index position in grid definition arrays assuming + # row-major flattening of arrays after numpy.meshgrid index_lat = ind / self.lon_size index_lon = ind % self.lon_size return d, index_lon.astype(np.int32), index_lat.astype(np.int32) diff --git a/tests/test_grid.py b/tests/test_grid.py index 08ad010..a5c84d1 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -137,6 +137,19 @@ def test_k_nearest_neighbor_list(self): assert gpi[1, 0] == 38430 assert gpi[1, 1] == 38429 + def test_nearest_neighbor_max_dist(self): + # test with maxdist higher than nearest point + gpi, dist = self.grid.find_nearest_gpi(14.3, 18.5, max_dist=100e3) + assert gpi == 25754 + assert len([dist]) == 1 + lon, lat = self.grid.gpi2lonlat(gpi) + assert lon == 14.5 + assert lat == 18.5 + + # test with maxdist lower than nearest point + gpi, dist = self.grid.find_nearest_gpi(14.3, 18.5, max_dist=10e3) + assert len(gpi) == 0 + assert len(dist) == 0 class TestCellGridNotGpiDirect(unittest.TestCase):