From fc917e05b86c711f4ed749ea7007cecf74147ff7 Mon Sep 17 00:00:00 2001 From: Sebastian Hahn Date: Wed, 10 Nov 2021 16:27:09 +0100 Subject: [PATCH] Add wgs84 conversion --- CHANGELOG.rst | 5 ++ environment.yml | 1 + src/fibgrid/construction.py | 80 ++++++++---------- ...n1650000.nc => fibgrid_sphere_n1650000.nc} | Bin ...d_n430000.nc => fibgrid_sphere_n430000.nc} | Bin ...n6600000.nc => fibgrid_sphere_n6600000.nc} | Bin src/fibgrid/realization.py | 41 ++++++++- tests/test_construction.py | 14 +-- tests/test_realization.py | 8 +- 9 files changed, 94 insertions(+), 55 deletions(-) rename src/fibgrid/files/{fibgrid_n1650000.nc => fibgrid_sphere_n1650000.nc} (100%) rename src/fibgrid/files/{fibgrid_n430000.nc => fibgrid_sphere_n430000.nc} (100%) rename src/fibgrid/files/{fibgrid_n6600000.nc => fibgrid_sphere_n6600000.nc} (100%) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 6a7bfa7..9beeef6 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -2,6 +2,11 @@ Changelog ========= +Version 0.0.2 +============= + +- Add Fibonacci grid realization for Geodatum WGS84 at 6.25, 12.5 and 25 km + Version 0.0.1 ============= diff --git a/environment.yml b/environment.yml index 53210f1..1756f26 100644 --- a/environment.yml +++ b/environment.yml @@ -5,6 +5,7 @@ channels: dependencies: - numpy - netCDF4 + - pyproj - numba - pip - pip: diff --git a/src/fibgrid/construction.py b/src/fibgrid/construction.py index d8209c5..d9f3f9c 100644 --- a/src/fibgrid/construction.py +++ b/src/fibgrid/construction.py @@ -32,7 +32,7 @@ import netCDF4 import numpy as np from numba import jit -from pygeogrids.grids import BasicGrid +from pyproj import CRS, Transformer @jit(nopython=True, cache=True) @@ -58,8 +58,8 @@ def compute_fib_grid(n): """ points = np.arange(-n, n+1) gpi = np.arange(points.size) - lat = np.empty(points.size) - lon = np.empty(points.size) + lat = np.empty(points.size, dtype=np.float64) + lon = np.empty(points.size, dtype=np.float64) phi = (1. + np.sqrt(5))/2. for i in points: @@ -73,52 +73,36 @@ def compute_fib_grid(n): return points, gpi, lon, lat -def compute_grid_stats(n, k=5, geodatum='WGS84'): +def compute_fib_grid_wgs84(n): """ - Compute grid statistics of k nearest neighbors. + Computation of Fibonacci lattice on a sphere and coordinated transformed + to WGS84 ellipsoid. Parameters ---------- n : int Number of grid points in the Fibonacci lattice. - k : int, optional - Nearest neighbor (default: 5). - geodatum : str, optional - Geodatum (default: 'WGS84') """ - points, gpi, lon, lat = compute_fib_grid(n) - - grid = BasicGrid(lon, lat, geodatum=geodatum) - nn, dist = grid.find_k_nearest_gpi(lon, lat, k=5) - - print('Min distance: {:10.4f}'.format(dist[:, 1:].min())) - print('Max distance: {:10.4f}'.format(dist[:, 1:].max())) - print('Mean distance: {:10.4f}'.format(dist[:, 1:].mean())) - print('Median distance: {:10.4f}'.format(np.median(dist[:, 1:]))) - print('Std distance: {:10.4f}'.format(dist[:, 1:].std())) - - for i in range(1, k): - print('\n') - print('------- Neighbor #{:} -------'.format(i)) - print('Min distance: {:10.4f}'.format(dist[:, i].min())) - print('Max distance: {:10.4f}'.format(dist[:, i].max())) - print('Mean distance: {:10.4f}'.format(dist[:, i].mean())) - print('Median distance: {:10.4f}'.format(np.median(dist[:, i]))) - print('Std distance: {:10.4f}'.format(dist[:, i].std())) - - print('\n') - lat_bands = range(-90, 90, 5) - for band in lat_bands: - subgrid = grid.subgrid_from_gpis(grid.get_bbox_grid_points( - latmin=band, latmax=band+5)) - nn, dist_subgrid = grid.find_k_nearest_gpi( - subgrid.arrlon, subgrid.arrlat, k=k) - print('Band [{} -- {}] deg: {:10.4f}'.format( - band, band+5, dist_subgrid[:, 1:].mean())) - - -def write_grid(filename, n, nc_fmt='NETCDF4_CLASSIC', nc_zlib=True, - nc_complevel=2): + crs_wgs84 = CRS.from_epsg(4326) + crs_sphere = CRS.from_proj4( + '+proj=lonlat +ellps=sphere +R=6370997 +towgs84=0,0,0') + + points, gpi, sphere_lon, sphere_lat = compute_fib_grid(n) + transformer = Transformer.from_crs(crs_sphere, crs_wgs84) + + wgs84_lon = np.zeros(sphere_lon.size, dtype=np.float64) + wgs84_lat = np.zeros(sphere_lat.size, dtype=np.float64) + + i = 0 + for lon, lat in zip(sphere_lon, sphere_lat): + wgs84_lat[i], wgs84_lon[i] = transformer.transform(lon, lat) + i = i + 1 + + return points, gpi, wgs84_lon, wgs84_lat + + +def write_grid(filename, n, nc_fmt='NETCDF4', nc_zlib=True, + nc_complevel=2, geodatum='WGS84'): """ Write grid file for Fibonacci lattice. @@ -137,7 +121,12 @@ def write_grid(filename, n, nc_fmt='NETCDF4_CLASSIC', nc_zlib=True, The optional keyword complevel is an integer between 1 and 9 describing the level of compression desired (default 2). """ - points, gpi, lon, lat = compute_fib_grid(n) + if geodatum == 'WGS84': + points, gpi, lon, lat = compute_fib_grid_wgs84(n) + elif geodatum == 'sphere': + points, gpi, lon, lat = compute_fib_grid(n) + else: + raise ValueError('Geodatum unknown') with netCDF4.Dataset(filename, 'w', format=nc_fmt) as fp: @@ -201,3 +190,8 @@ def read_grid(filename, variables=['gpi', 'lon', 'lat']): data[var] = fp.variables[var][:] return data + + +if __name__ == '__main__': + n = 6600000 + compute_fib_grid_wgs84(n) diff --git a/src/fibgrid/files/fibgrid_n1650000.nc b/src/fibgrid/files/fibgrid_sphere_n1650000.nc similarity index 100% rename from src/fibgrid/files/fibgrid_n1650000.nc rename to src/fibgrid/files/fibgrid_sphere_n1650000.nc diff --git a/src/fibgrid/files/fibgrid_n430000.nc b/src/fibgrid/files/fibgrid_sphere_n430000.nc similarity index 100% rename from src/fibgrid/files/fibgrid_n430000.nc rename to src/fibgrid/files/fibgrid_sphere_n430000.nc diff --git a/src/fibgrid/files/fibgrid_n6600000.nc b/src/fibgrid/files/fibgrid_sphere_n6600000.nc similarity index 100% rename from src/fibgrid/files/fibgrid_n6600000.nc rename to src/fibgrid/files/fibgrid_sphere_n6600000.nc diff --git a/src/fibgrid/realization.py b/src/fibgrid/realization.py index b68aaec..1d25ae5 100644 --- a/src/fibgrid/realization.py +++ b/src/fibgrid/realization.py @@ -36,7 +36,7 @@ from pygeogrids.grids import CellGrid -def read_grid_file(n): +def read_grid_file(n, geodatum='WGS84'): """ Read pre-computed grid files. @@ -59,8 +59,9 @@ def read_grid_file(n): metadata : dict Metadata information of the grid. """ - filename = os.path.join(os.path.dirname(__file__), - 'files', 'fibgrid_n{}.nc'.format(n)) + filename = os.path.join( + os.path.dirname(__file__), 'files', + 'fibgrid_{}_n{}.nc'.format(geodatum.lower(), n)) metadata_fields = ['land_frac_fw', 'land_frac_hw', 'land_mask_hw', 'land_mask_fw'] @@ -104,5 +105,37 @@ def __init__(self, res, geodatum='WGS84'): else: raise ValueError('Resolution unknown') - lon, lat, cell, gpi, self.metadata = read_grid_file(n) + lon, lat, cell, gpi, self.metadata = read_grid_file( + n, geodatum=geodatum) + super().__init__(lon, lat, cell, gpi, geodatum=geodatum) + + +class FibLandGrid(CellGrid): + + """ + Fibonacci grid with active points over land defined by land fraction. + """ + + def __init__(self, res, geodatum='WGS84'): + """ + Initialize FibGrid. + + Parameters + ---------- + res : int + Sampling. + geodatum : str, optional + Geodatum (default: 'WGS84') + """ + if res == 6.25: + n = 6600000 + elif res == 12.5: + n = 1650000 + elif res == 25: + n = 430000 + else: + raise ValueError('Resolution unknown') + + lon, lat, cell, gpi, self.metadata = read_grid_file( + n, geodatum=geodatum) super().__init__(lon, lat, cell, gpi, geodatum=geodatum) diff --git a/tests/test_construction.py b/tests/test_construction.py index 1833fb5..2646598 100644 --- a/tests/test_construction.py +++ b/tests/test_construction.py @@ -46,6 +46,7 @@ def setUp(self): """ self.test_path = mkdtemp() self.n = [6600000, 1650000, 430000] + self.geodatum = ['sphere', 'WGS84'] def test_fibgrid(self): """ @@ -60,11 +61,14 @@ def test_read_write_fibgrid(self): """ Test read/write Fibonacci grid. """ - for n in self.n: - filename = os.path.join(self.test_path, 'fibgrid_{}.NC'.format(n)) - write_grid(filename, n) - data = read_grid(filename) - np.testing.assert_equal(data['gpi'], np.arange(n*2+1)) + for geodatum in self.geodatum: + for n in self.n: + filename = os.path.join( + self.test_path, 'fibgrid_{}_n{}.nc'.format( + geodatum.lower(), n)) + write_grid(filename, n) + data = read_grid(filename) + np.testing.assert_equal(data['gpi'], np.arange(n*2+1)) if __name__ == '__main__': diff --git a/tests/test_realization.py b/tests/test_realization.py index 289dcf9..57023be 100644 --- a/tests/test_realization.py +++ b/tests/test_realization.py @@ -42,14 +42,16 @@ def setUp(self): """ self.res = [6.25, 12.5, 25] self.n = [6600000, 1650000, 430000] + self.geodatum = ['sphere', 'WGS84'] def test_grid_files(self): """ Test read grid files. """ - for res, n in zip(self.res, self.n): - fb = FibGrid(res) - assert fb.gpis.size == 2*n+1 + for geodatum in self.geodatum: + for res, n in zip(self.res, self.n): + fb = FibGrid(res, geodatum=geodatum) + assert fb.gpis.size == 2*n+1 if __name__ == '__main__':