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

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

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