I am using R to loop through thousands of NetCDF files (downscaled climate data), and while I am able to do it using R, the code is too slow. I am looking for a way to speed it up and to maintain functionality.

What I need is a program or platform where I can (1) import netcdf data, (2) import a polygon layer, (3) get the mean for the polygon layer, and (4) export output to csv file. The following is a general code workflow:

shape <- readOGR(dsn,"shape") ## read in shapefile
files <- dir("in directory", recursive=TRUE, full.names=TRUE, pattern="\\.nc$") ### get file paths to netCDF files
out <-NULL
for (i in 1:length(files)){
brick.tmp <- files[i]
    for (j in 1:nlayers(brick.tmp)){
    val = extract(brick.tmp[[j]],shape, fun = mean)
    out <- rbind(out,val)

I would prefer to stay in the R environment, but I'm not sure how to speed up the code.

I have attempted to use several programs and have looked through the list here: https://www.unidata.ucar.edu/software/netcdf/docs/software.html

However, I thought I would leverage everyone's collective knowledge and see if I can narrow down the potential choices.

  • Which dimensions does your netcdf file has? I would suspect X,Y,T? Do you want to get the mean over all cells within the polygon and over time? – ulrich Sep 3 '15 at 13:59
  • 1
    Before jumping out to a completely different language perhaps you should explore ways to speed up your R code. However, since you do not provide any details on what you have tried in R we can give no insight as to ways to potentially speed up you code. As the saying goes "there is more than one way to skin a cat" which, especially applies to R. – Jeffrey Evans Sep 3 '15 at 15:39
  • Thanks you two. it is indeed X,Y,T. @Jeffrey Evans, you are correct and I will edit my question to take into account your comments. – user44796 Sep 3 '15 at 15:51
  • 1
    I have found that for zonal like functions it is best to actually mask your raster(s) to the polygon and then just pull the statistic out of the results by coercing to a vector or matrix. You will be surprised how much faster this is than extract. I have run some large problems (10,000 polygons with 0.97m raster). When you mask a stack it does all rasters at once so, there is no real performance ding. – Jeffrey Evans Sep 3 '15 at 16:08
  • @JeffreyEvans Are you talking about looping over the polygons to avoid using extract? If you could provide a few example lines of code, that would be great. Thanks. – user44796 Sep 3 '15 at 18:43

NetCDF is incredibly general and writing slow code is easy. I routinely deal with tens of thousands of NetCDF files in R, using some combination of packages raster, ncdf, ncdf4, RNetCDF, or rgdal. The key is to leverage the cell index tools in raster so that the "cell-in-polygon" test occurs only once, then you can apply the extraction across all files. Many of the high-level tools do this for you so it can be very easy.

A pseudo workflow looks like:

## build a sensible file database (more than this, amazing how rarely it is done)
fs <- list.files("nctree", etc. etc.)  
poly <- readOGR("folder/of/polys", "shapefilename")
polyMatrixVals <- extract(stack(fs), poly, fun = mean)

Lots of assumptions are loaded in there, but if it doesn't quite fit there are ways to get the same effect using the component tools in raster. Whether R is the right tool depends on the details of your situation, and that's true for any software choice.

  • @mdsummer You're recommendation reduced my runtime on a subset of the data from 71 minutes (by looping through the brick layers) to 28 minutes. That is a step in the right direction. – user44796 Sep 3 '15 at 18:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.