How to extract value from a raster based on lat and lon?

Question

I have several netcdf files that I can read using the netCDF4 library in python. The data is gridded in 4320 rows and 2160 columns.

I also have a list of lat/lon (in decimal degrees) and I am trying to find the value of the gridcell in which the coordinates lie.

Currently, I am converting the coordinates into gridcell boundaries by a simple linear transformation and querying the netcdf file using these new points. Is this approach valid?

# prelim
from netCDF4 import Dataset
import numpy as np

# coordinate data
lons = [40.822651, 40.69685, 40.750038, 40.798931, 40.74612]
lats = [19.83832, 20.007561, 19.97426, 19.86334, 19.843889]

# find boundaries of gridcells corresponding to lat/lon
lons_grid = np.ceil(12*(lons + 180))
lats_grid = np.ceil(12*(lats + 90))

# find value in gridcell for given lat/lon
fnc = Dataset(filename.nc, 'r')
lat = fnc.variables['latitude'][:]
lon = fnc.variables['longitude'][:]
variable = fnc.variables['variable'][:]
print variable[lons_grid, lats_grid]

If not, is there a way to intersect point data with a raster in python without using arcpy? (My netcdf files are very large and cannot be read with arcpy.)

Answer 1

To get the appropriate cell coordinates from your latitude and longitude you need to know the coverage of the netCDF file (for example, the coordinate of the upper left cell and the size of the cells in the x and y directions). @Luke 's solution below works for a srtaight x/y raster with no cell rotation. If you don't have that, you could look at the affine library to transform latitude and longitude to cell coordinates. Given GDAL GeoTransform in the form:

1. Upper left cell x coordinate
2. x cell size in projected units
3. x directional rotation (usually 0)
4. Upper left cell y coordinate
5. y directional rotation (usually 0)
6. Negative y cell size in projected units

So for example the -237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0 would translate as an upper left cell of the raster with coordinates -237481.5, 237536.4, and a 425.0 unit square with no rotation.

Using the affine library you can transform this into a transformation object like so:

from affine import Affine
aff = Affine.from_gdal(-237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0)

Which can transform cell coordinates to projected coordinates:

xs = np.array([0, 1, 2])
ys = np.array([3, 4, 5])

x_coords, y_coords = aff * (xs, ys)

Or (in your case) from coordinates back to cell coordinates:

xs, ys = ~aff * (np.array(lons), np.array(lats))

These values are floats, so you'll need to transform them into integers to get cell coordinates you can use.

xs = np.round(xs).astype(np.int)
ys = np.round(ys).astype(np.int)

You can then use these as indexes to your netCDF4 array (use the latest version of netCDF4 - 1.2.1 as this means you don't need to sort or remove duplicates).

variable = fnc.variables['variable'][xs, ys] # be careful of your index order

This will return a square array due to the slightly different netCDF indexing, but you can get the actual values that you're after as the diagnoal:

values = variable.diagonal()

Answer 2

If you know the pixel/cell size and minimum lon/lat of the grid, it's a simple calculation.

px = (lon -  minlon) / xcellsize
py = (lat -  minlat) / ycellsize

You can also use min/max lon/lat and the number of cols/rows:

xcellsize =  (maxlon -  minlon) / ncols
ycellsize =  (maxlat -  minlat) / nrows
#then use previous calculation to get pixel coordinates