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:
- Upper left cell
x coordinate
x cell size in projected unitsx directional rotation (usually 0)- Upper left cell
y coordinate
y directional rotation (usually 0)- 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()
Here's a way, first I'll create a fake data set.
library(raster)
r <- raster(matrix(1:30, 5, 6))
## this is the full data set in data frame form
dfull <- as.data.frame(r, xy = TRUE)
## this is the partial data set, only the points with a valid value
## row-order doesn't matter, but we keep it for illustration
set.seed(10)
dpart <- dfull[sort(sample(seq_len(nrow(dfull)), 22)), ]
Now we need raster's cell-abstraction tools. Here we can treat r like an raw specification of the original raster, and in fact create it from scratch if needed. But, we have it so we use it.
rspec <- raster(r) ## this drops the data, keeps the structure
This fills the data with missing values , because raster doesn't truly have "sparse forms", they are either empty or full and we cannot put values piece-wise into an empty raster, it's either all or nothing until it's not empty.
(Note that sparse forms are supported completely by this approach, but via a level of abstraction that is the responsibility of the user)
rspec[] <- NA_real_
Now we need an index into the "structure of the raster" for our points.
## these names were nominated above, and might be different for a different
## input
cells <- cellFromXY(rspec, as.matrix(dpart[, c("x", "y")])
Now, put the values in the data frame into the otherwise "filled with missing" raster.
rspec[cells] <- dpart$layer
All this is illustrated in full here: http://rpubs.com/cyclemumner/294656
It's a very powerful approach, but it's not widely understood and it's easy to get it wrong, so do use with caution and take time to practice and understand it.
Happy to help if it doesn't make sense.
Best Answer
You can do
But your question suggests that you would have been better of accessing the ncdf file as a RasterLayer object (and avoid creating matrices and keeping track of coordinates)
etc.