Problem with raster package and nc4 file

You don't need any hacky workarounds, this is a respectable NetCDF file. I
have an account on that system and found a file called "tmin_1999.nc4" -
this is large, ~ 1.5Gb.   Extracting pixel values with a polygon object is
pretty straightforward, but this is quite a large data set so

- test carefully on just a single layer
- check that the automatic transformation of your polygons to the grid by
extract() is sensible (or do it yourself)
- consider subsetting the grid if you don't need all of its extent in
space/time
- extract will be efficient for multilayer grids, but you might consider
use cellnumbers = TRUE if you want to save the index for later extractions

HTH

Cheers, Mike.

## this requires registration at
http://daac.ornl.gov/DAYMET/guides/Daymet_mosaics.html
##f <- "
ftp://ftp.daac.ornl.gov/data/daymet/Daymet_mosaics/data/tmin_1999.nc4"
##download.file(f, basename(f), mode = "wb")
## (could also use OpenDAP to read this if you have rgdal built with that)

## all the information is here:
##http://daac.ornl.gov/DAYMET/guides/Daymet_mosaics.html

## raster already knows:
r <- stack( basename(f))
class       : RasterStack
dimensions  : 4823, 5268, 25407564, 365  (nrow, ncol, ncell, nlayers)
resolution  : 1000, 1000  (x, y)
extent      : -2015000, 3253000, -3037000, 1786000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +lat_1=25
+lat_2=60 +ellps=WGS84
names       : X1999.01.01, X1999.01.02, X1999.01.03, X1999.01.04,
X1999.01.05, X1999.01.06, X1999.01.07, X1999.01.08, X1999.01.09,
X1999.01.10, X1999.01.11, X1999.01.12, X1999.01.13, X1999.01.14,
X1999.01.15, ...

## check with a map
library(rworldmap)
data(countriesLow)
library(rgdal)
namerica <- subset(countriesLow, SOVEREIGNT %in% c("United States of
America", "Mexico", "Canada"))
wm <- spTransform(namerica, CRS(projection(r)))

plot(r[[1]])
plot(wm, add = TRUE)  ## looks fine to me . . .

## this code untested obviously
## read in your shapefile
# poly <- shapefile("poly.shp")
# polyvals <- extract(r, poly, fun = mean, na.rm = TRUE) ## use r[[1]] to
test on a single layer and see notes above

## dumb polygon extraction test - this will take some time to run, so try
small tests first
extract(r[[1]], wm, fun = mean, na.rm = TRUE)

Problem with raster package and nc4 file

Thread (6 messages)