I'm trying to display mean annual precipitation (MAP) data from the CRU dataset (http://data.ceda.ac.uk/badc/cru/data/cru_ts/cru_ts_4.01/data/pre/) across a specific region.

The data are monthly precipitation totals. To get MAP, I need to first calculate annual sum for each year (1901-2016), then average across all years. Ideally, I would also subset this into years 1950-2013. The CRU data come in netCDF format, which I'm unfamiliar with. Do I need to convert the netCDF to a dataframe, make the MAP calculations, then convert back to netCDF? My ultimate goal is to then export the MAP data in gridded format as a raster to ArcMap.

Below is the R code I've been using for overlaying the overall average CRU data on my polygon. However, the values created by prec_rasta_mean below are not MAP.


## Load shapefile to R
shp <- shapefile("GridAK_IU2016_FINAL_PolyAA83_SCSE")

## Plot shape to make sure

## Retrieve data within shape file across timeRange
prec_rasta <- cruts2raster("cru_ts4.01.1901.2016.pre.dat.nc", 
timeRange=c("1950-01-01", "2013-12-31"), shp, type='brick')

## Average data across all months in timeRange
prec_rasta_mean <- calc(prec_rasta, mean)

##  Write new raster layer to wd
writeRaster(prec_rasta_mean, 'cru_prec_brick_mean.img')

So, I'm looking for a bit of help getting making the calculations in netCDF, then I think the code above should take care of the rest.

The polygon that I'm using was created in ArcMap and the shapefile is available here: https://drive.google.com/file/d/0B_eqTercwIH2d01tekdaY2tSbEU/view?usp=sharing

Bit of a novice playing with spatial data in R (and ArcMap).


1 Answer 1


I came up with a solution (with help from colleagues).

The following code gets MAP for each grid cell, averaged over the timerange of interest (1950-2013). Basically I created a vector that could be used to index each year, summed across years, then took the average of the whole timerange.

# create vector to serve as index
num_years <- rep(1:113, each=12)

# sum annual precipitation across months for each year, indexed using 
# num_years
rasta_sum <- stackApply(prec_rasta, indices=num_years, fun=sum)

# average annual sum for each cell across years
rasta_mean <- calc(rasta_sum, mean, na.rm=TRUE)

In the end, I didn't need to manipulate the netCDF file at all. Just converted to raster using cruts2raster function, then performed calculations with raster.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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