Basics

import pygeogrids.grids as grids
import numpy as np

Let’s create a simple regular 10x10 degree grid with grid points at the center of each 10x10 degree cell.

First by hand to understand what is going on underneath

# create the longitudes
lons = np.arange(-180 + 5, 180, 10)
print(lons)
lats = np.arange(90 - 5, -90, -10)
print(lats)

[-175 -165 -155 -145 -135 -125 -115 -105  -95  -85  -75  -65  -55  -45
  -35  -25  -15   -5    5   15   25   35   45   55   65   75   85   95
  105  115  125  135  145  155  165  175]
[ 85  75  65  55  45  35  25  15   5  -5 -15 -25 -35 -45 -55 -65 -75 -85]

These are just the dimensions or we can also call them the “sides” of the array that defines all the gridpoints.

# create all the grid points by using the numpy.meshgrid function
longrid, latgrid = np.meshgrid(lons, lats)

now we can create a BasicGrid. We can also define the shape of the grid. The first part of the shape must be in longitude direction.

manualgrid = grids.BasicGrid(longrid.flatten(), latgrid.flatten(), shape=(18, 36))

# Each point of the grid automatically got a grid point number
gpis, gridlons, gridlats = manualgrid.get_grid_points()
print(gpis[:10], gridlons[:10], gridlats[:10])

[0 1 2 3 4 5 6 7 8 9] [-175 -165 -155 -145 -135 -125 -115 -105  -95  -85] [85 85 85 85 85 85 85 85 85 85]

The grid point indices or numbers are useful when creating lookup tables between grids.

We can now use the manualgrid instance to find the nearest gpi to any longitude and latitude

ngpi, distance = manualgrid.find_nearest_gpi(15.84, 28.76)
print(ngpi, distance)
# convert the gpi to longitude and latitude
print(manualgrid.gpi2lonlat(ngpi))

235 424808.5131778209
(15, 25)

The same grid can also be created by a method for creating regular grids

autogrid = grids.genreg_grid(10, 10)
autogrid == manualgrid

True

If your grid has a 2D shape like the ones we just created then you can also get the row and the column of a grid point. This can be useful if you know that you have data stored on a specific grid and you want to read the data from a grid point.

row, col = autogrid.gpi2rowcol(ngpi)
print(row, col)

6 19

Iteration over gridpoints

for i, (gpi, lon, lat) in enumerate(autogrid.grid_points()):
    print(gpi, lon, lat)
    if i==10: # this is just to keep the example output short
        break

0 -175.0 85.0
1 -165.0 85.0
2 -155.0 85.0
3 -145.0 85.0
4 -135.0 85.0
5 -125.0 85.0
6 -115.0 85.0
7 -105.0 85.0
8 -95.0 85.0
9 -85.0 85.0
10 -75.0 85.0

Calculation of lookup tables

If you have a two grids and you know that you want to get the nearest neighbors for all of its grid points in the second grid you can calculate a lookup table once and reuse it later.

# lets generate a second grid with 10 random points on the Earth surface.

randlat = np.random.random(10) * 180 - 90
randlon = np.random.random(10) * 360 - 180
print(randlat)
print(randlon)
# This grid has no meaningful 2D shape so none is given
randgrid = grids.BasicGrid(randlon, randlat)

[-30.89674861  59.08649779 -49.27545779 -69.52856569  -1.39290501
 -10.88163146  53.67514799  60.14998259  55.15218786 -29.79836856]
[ 125.68580631 -169.379169    105.7270734    46.42576672  -28.28325787
   47.49761201  173.67561423   86.59161319 -178.39301312  115.68242529]

Now lets calculate a lookup table to the regular 10x10° grid we created earlier

lut = randgrid.calc_lut(autogrid)
print(lut)

[462 109 496 562 339 382 143  98 108 425]

The lookup table contains the grid point indices of the other grid, autogrid in this case.

lut_lons, lut_lats = autogrid.gpi2lonlat(lut)
print(lut_lats)
print(lut_lons)

[-35.  55. -45. -65.  -5. -15.  55.  65.  55. -25.]
[ 125. -165.  105.   45.  -25.   45.  175.   85. -175.  115.]

Storing and loading grids

Grids can be stored to disk as CF compliant netCDF files

import pygeogrids.netcdf as nc
nc.save_grid('example.nc', randgrid)

loadedgrid = nc.load_grid('example.nc')

loadedgrid

<pygeogrids.grids.BasicGrid at 0x7f4110fe7700>

randgrid

<pygeogrids.grids.BasicGrid at 0x7f41128d92b0>

Define geodetic datum for grid

grid_WGS84 = grids.BasicGrid(randlon, randlat, geodatum='WGS84')

grid_GRS80 = grids.BasicGrid(randlon, randlat, geodatum='GRS80')

grid_WGS84.geodatum.geod.a

6378137.0

grid_GRS80.geodatum.geod.a

6378137.0

grid_WGS84.kdTree.geodatum.geod.sphere

False