# Using Proj.4/GDAL to convert lat,lon to i,j

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
• 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
``````