[R-sig-Geo] get ncdf values at point locations - efficient method?

[Robert J. Hijmans]

Noah,
There is a third way:

library(raster)
b <- brick("test.nc")

(this does not load the cell values into memory). Now there are
several ways to extract subsets. See ?getValues, ?crop, ....
Robert



On Sun, Nov 17, 2013 at 9:40 PM,  <ndcNS at hampshire.edu> wrote:
> Hi,
>
> I'm trying to extract data from point locations in NetCDF files efficiently.
> I can do it for small data, but I am moving to scenarios with huge grids and
> 500,000 points.  For small data, I do one of two things:
>
> 1) Read the entire grid into memory and convert it to a raster stack, then
> extract cell values at my coordinates.
> 2) Loop through each point, reading in just the single value from the ncdf.
>
> The problem with #1 is that it is too slow and large to fit all those
> extraneous grid cells into memory.  The problem with #2 is that looping
> through 500,000 points is too slow.  What I want to do is something like #2,
> but without a loop - a single vectorizable call.  Any suggestions?
>
> Here's a simplified version of what my #2 code looks like now:
>
> opened.ncdf <- open.ncdf(test.nc)
> lons <- get.var.ncdf(opened.ncdf, 'longitude')
> lats <- get.var.ncdf(opened.ncdf, 'latitude')
> point.coords <- read.csv('lonlat.csv')
> point.coords[,1] <- sapply(point.coords[,1],function(x)which.min(abs(x -
> lons)))
> point.coords[,2] <- sapply(point.coords[,2],function(x)which.min(abs(x -
> lats)))
> layer.num <- 3 # this is time, I have seperate ways for handling this, let's
> just call it timestep 3 for now
> output <- vector('numeric',length(dim(point.coords)[1]))
>
> for(i in 1:(dim(point.coords)[1])){
>         output[i] <- get.var.ncdf(opened.ncdf,
> start=c(point.coords[i,],layer.num), count=c(1,1,1))
> }
>
> Is there a way to vectorize that slow loop???
>
> Thanks!!!
> -Noah
>
> ---
> Noah Charney, PhD
> Visiting Assistant Professor of Ecology
> School of Natural Science
> Hampshire College
> 304 Cole Science Center
> Amherst, MA 01002
> 413-559-5775
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo