[R-sig-Geo] ncdf 4-dim file to raster brick - level error??

Gross, John john_gross at nps.gov
Sun Oct 4 21:48:16 CEST 2015


I'm trying to read a 4-dimensional ncdf file (dims: lat lon time
projection) into a raster brick and can't get it to read anything
other than the first "level" (i.e. projection, the 4th dimension).

Here's the code, output, and sessionInfo.  This is reminiscent of a
problem a couple years ago with earlier versions of both ncdf4 and
raster that was caused by the ncdf file not closing between calls  But
I'm not sure how to determine it that's happening again.

You can download the data file (10KB - it's only 3x3x59x4) from:
https://www.dropbox.com/s/5jh5n8etpk1y1hu/Extraction_tasmin.nc?dl=0

#===  Code ==========
library(ncdf4)
library(raster)

nc <- nc_open("Extraction_tasmin.nc")    # this is a 3x3 text grid
nc    # dims are lat long time and projection

xy <- data.frame(x = min(nc$dim$lon$vals) + 1/32, y =
max(nc$dim$lat$vals) - 1/32) #get top left point not on grid boundary
coordinates(xy) <- ~x+y

for(i in 1:5){
  b1 <- brick("Extraction_tasmin.nc", lvar=4, level=i)
  pts <- t(extract(b1, xy))
  print(c("pts[1:5} -> ",pts[1:5])); flush.console()
}
#====  Output  =======

File Extraction_tasmin.nc (NC_FORMAT_CLASSIC):

     1 variables (excluding dimension variables):
        float tasmin[longitude,latitude,time,projection]
            typeConversion_op_ncl: double converted to float
            _FillValue: 1.00000002004088e+20

     4 dimensions:
        latitude  Size:3
            valid_max: 52.875
            long_name: Latitude
            valid_min: 25.125
            units: degrees_north
            axis: Y
        longitude  Size:3
            long_name: Longitude
            axis: X
            units: degrees_east
            modulo: 360
            topology: circular
        time  Size:59
            calendar: standard
            units: days since 1950-01-01 00:00:00
            standard_name: time
            long_name: time
            axis: T
        projection  Size:4   *** is unlimited ***

    2 global attributes:
        NCO: 4.0.8
        Projections: access1-0.1.rcp85, bcc-csm1-1.1.rcp85,
canesm2.1.rcp85, ccsm4.1.rcp85,
>

Output:
[1] "pts[1:5} -> "      "-5.08642387390137" "-3.23856282234192"
"-5.46992635726929" "-3.83646535873413" "-3.00846982002258"
[1] "pts[1:5} -> "      "-5.08642387390137" "-3.23856282234192"
"-5.46992635726929" "-3.83646535873413" "-3.00846982002258"
[1] "pts[1:5} -> "      "-5.08642387390137" "-3.23856282234192"
"-5.46992635726929" "-3.83646535873413" "-3.00846982002258"
[1] "pts[1:5} -> "      "-5.08642387390137" "-3.23856282234192"
"-5.46992635726929" "-3.83646535873413" "-3.00846982002258"
[1] "pts[1:5} -> "      "-5.08642387390137" "-3.23856282234192"
"-5.46992635726929" "-3.83646535873413" "-3.00846982002258"
Warning message:
In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
  level set to: 4
# ==== sessionInfo =========
R version 3.2.2 (2015-08-14)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] ncdf4_1.13    raster_2.4-18 sp_1.1-1

loaded via a namespace (and not attached):
[1] rgdal_1.0-4     tools_3.2.2     Rcpp_0.12.1     grid_3.2.2
lattice_0.20-33

-- 
John Gross, PhD
NPS, Ft Collins, CO



More information about the R-sig-Geo mailing list