I'm working with a data set in GRIB2 format and need to determine the grid cell index (i, j) corresponding to a given lat,lon. I have the projection definition in the form of projparams as follows:

{'a': 6371229,
 'b': 6371229,
 'lat_0': 25.0,
 'lat_1': 25.0,
 'lat_2': 25.0,
 'lon_0': 265.0,
 'proj': 'lcc'}

which gives me the following srs:+units=m +a=6371229 +b=6371229 +lon_0=265.0 +proj=lcc +lat_2=25.0 +lat_1=25.0 +lat_0=25.0

And I know dx, dy, number of grid points in each direction etc.

So I was thinking of using Proj.4 (because it has JavaScript bindings) (and/or GDAL) to convert it to i,j. I think I have read somewhere that one could do it this way by defining some kind of equidistant metric map projection that would use i, j as x, y coordinates and do a simple coordinate transformation, but I'm struggling with how to actually define such a projection in proj.4 terms.

  • Can't you just use gdallocationinfo directly? – BradHards Dec 10 '15 at 9:00
  • @BradHards Maybe I can. A Proj.4 way seems more lightweight for this task. GDAL is a much bigger library. – exfizik Dec 11 '15 at 7:40
  • 1
    I just put that in because your question (including the original tags and the title) mentioned gdal. – BradHards Dec 11 '15 at 8:30
  • @BradHards, Yes, you're right. I wasn't sure which one I'd use at first. – exfizik Dec 11 '15 at 15:31

Here is how wgrib2 utility does it. The algorithm:

  1. Initialize 2 projections: an x,y one called pj_grid and a lat,lon one called pj_latlon
  2. Given lat,lon of the first grid cell, find its x,y coords (x_0, y_0)
  3. Given lat,lon of a point of interest, find its x,y coords (x,y)
  4. Given the grid spacings dx, dy, find the grid indices of the point of interest i,j as:

    i = (x-x_0)/dx + 1

    j = (y-y_0)/dy + 1

Note: the indices are 1 based. Below is a sample implementation in Python using pyproj. The grid used is NAM218 grid

import pyproj
# in my case projparams is obtained using pygrib
projparams = {
    'a': 6371229,
    'b': 6371229,
    'lat_0': 25.0,
    'lat_1': 25.0,
    'lat_2': 25.0,
    'lon_0': 265.0,
    'proj': 'lcc'

pj_grid = pyproj.Proj(projparams=projparams)
print pj_grid.srs 
# +units=m +a=6371229 +b=6371229 +lon_0=265.0 +proj=lcc +lat_2=25.0 +lat_1=25.0 +lat_0=25.0
pj_latlon = pj_grid.to_latlong()
print pj_latlon.srs
# +proj=latlong +a=6371229 +b=6371229

# grid spacing in m
dx = 12191
dy = 12191

# first grid point, lat/lon in degrees
lat1 = 12.19
lon1 = 226.541

x_0, y_0 = pyproj.transform(pj_latlon, pj_grid, lon1, lat1, radians=False)
print x_0, y_0    # -4226106.99692 -832698.261018

lat = 51
lon = -114

x, y = pyproj.transform(pj_latlon, pj_grid, lon, lat, radians=False)
# these are 1 based indices
i = (x - x_0) / dx + 1
j = (y - y_0) / dy + 1
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.