1

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.

library(cruts)
library(raster)
library(rgdal)

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

## Plot shape to make sure
plot(shp)

## 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).

  • Welcome to GIS SE! As a new user be sure to take the Tour to learn about our focussed Q&A format. I've removed your tag for ArcMap because your current question seems to be about where you are stuck with R. – PolyGeo Sep 30 '17 at 1:19
  • did you try rowMeans()? – Elio Diaz Oct 4 '17 at 19:04
0

Just returning to this project and 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. Anyways, hope this helps.

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.