Add wgs84 conversion

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
 =============
 

environment.yml

@@ -5,6 +5,7 @@ channels:
 dependencies:
   - numpy
   - netCDF4
+  - pyproj
   - numba
   - pip
   - pip:

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)

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)

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__':

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__':