2

I need to build a table that has monthly precipitation and temperature means for basically every country I can (from 2001-2014) using the NCEP/NCAR Reanalysis Monthly Means Dataset.

I have been able to load the temperature and precipitation monthly means data into both ArcGIS and R, along with a shapefile for the country borders from thematicmapping.org, but I have been completely unable to extract data for each month with either method.

In ArcGIS I have both datasets loaded and I can use the slider to see different months, as well as use the Zonal Statistics as Table tool, but I can't figure out how to get that to output every month. In R I have been able to load the files in but not get much deeper.

The NCEP/NCAR data is netcdf4 and can be found at: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.surface.html (air.mon.mean.nc and pr_wtr.eatm.mon.mean.nc).

The shapefile can be found at: http://thematicmapping.org/downloads/world_borders.php

  • It might be helpful to add a link to the data you are using (particularly the climate data) so we know what format it's in. I'm guessing it's a netcdf. It sounds like you're doing things correctly, you just need to go to the next step. Each month is going to be a different band of data, and zonal stats runs on one band at a time. So you'll need to make a script or model that iterates through each band, runs the zonal stats, and then appends that information in a new column/field to a table. – Chris W Apr 3 '15 at 19:20
  • Thanks for taking the time to write that up. I've added the exact data I'm trying to use. I understand the basic principle of what I'm trying to accomplish but I have no idea how to accomplish it in the model builder. I'm sure this is impossibly simple but do you happen to know of any model examples I could work from? – Spencer Dorsey Apr 3 '15 at 19:35
  • Unfortunately I'm not very familiar netcdfs and so I don't know the tools and methods to extract or address their values (let alone which ones are best or if there is a better way to work with them than the raster band concept I mention above). Hopefully mdsumner's answer will at least get you headed in the right direction. – Chris W Apr 3 '15 at 23:06
3

We have one of these files on hand:

## dp is the root to our local data repository
f <- file.path(dp, "data", "ftp.cdc.noaa.gov/Datasets/ncep.reanalysis2.derived/pressure/air.mon.mean.nc")
library(raster)
b <- brick(f, level = 1)
b

On my system that uses the "ncdf4" package, but it could also use "ncdf".

The global Pacific-view extent looks close, but I would reset it after checking carefully in the file:

extent : -1.25, 358.75, -91.25, 91.25 (xmin, xmax, ymin, ymax)

library(ncdf4)
nc <- nc_open(f)
## print(nc) to see details 
lon <- ncvar_get(nc, "lon")
range(lon)
##[1]   0.0 357.5

So, the longitudes that start at 0 and end somewhere near 0 need to be fixed, and so do the latitudes. I think it's pretty uncontroversial to assume they mean to fill the whole world (raster makes the reasonable assumption that each coordinate belongs to the centre, but clearly that's not what these are. It's typical for these incomplete specifications to exist, you are expected to fix it yourself. NetCDF has no standard concept of an affine transformation afaics so it's stuck with these redundant and often broken 1D coordinate arrays for regular grids.

extent(b) <- extent(0, 360, -90, 90)

The function rotate() will put the object into Atlantic view.

You can read other levels ("air" is a 4D variable with 432 time steps atm and 17 height levels) by setting that in the brick() argument.

I tend to always work this way:

  • use raster to see what happens

  • read the metadata with some netcdf package (ncdf, ncdf4, RNetCDF, or possibly GDAL or ncdump on the command line, occasionally Manifold GIS)

  • use those investigations to hone the use of raster.

You can't really use raster for non-regular coordinates bound to arrays, unless you really know what you're doing (but this air temp data is fine for XYT). See the available levels with

ncvar_get(nc, "level")
#[1] 1000  925  850  700  600  500  400  300  250  200  150  100   70   50   30   20   10  

10 is level = 1 for raster, etc.

Using raster you can get to individual months like this:

b[[10]]

subset(b, 10:13)

but now we get into "go learn the raster package" territory.

I use the "CountriesLow" data set in "rworldmap" to get mean country temperatures:

asub <- which(as.Date(getZ(b)) >= as.Date("2001-01-01"))

library(rworldmap)
data(countriesLow)
## extract the first year's data
imonth <- 1
d <- data.frame(Name = countriesLow$SOVEREIGNT, meanairtemp = extract(rotate(b[[asub[imonth]]]), countriesLow, fun = mean, na.rm = TRUE),
            date = as.Date(getZ(b)[asub[imonth]]), stringsAsFactors = FALSE)

I've assumed you would build a long-form table and just append each year to this, but you might want a column per month? One of the values comes out as NA but I think that can be controlled by small/sampling options with extract (your polygon layer may differ).

Don't forget to rotate() to get the Pacific view to match the Atlantic polygons, and that the data is in Kelvin.

Use rgdal::readOGR to get your shapefile or whatever in. extract() will automatically reproject from different projections, but that won't always work or apply (certainly not for the Pacific/Atlantic conventions).

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.