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.)