[R-sig-Geo] raster working with netCDF files that have indexed geography

[Anthony Fischbach]

In working with Global Circulation Model datasets from the Earth System Grid
Federation (http://pcmdi9.llnl.gov/esgf-web-fe/), and am struggling to
efficiently pull in the geography of the netCDF files.
My first efforts are based on the sea ice concentration (variable = sic)
projected by the NCAR rcp8.5 through the end of this century.

Optimally we should be able to follow the eloquent code available in the
raster package to read in the files.
## Create a raster brick from the netCDF 
b<-brick(x='filesic_OImon_CCSM4_rcp85_r1i1p1_200601-210012.nc',
varname='sic') 
##
This quickly reads in the netCDF file, however, it fails to read in the true
geography.
Even as it reports the projection as geographic (+proj=longlat
+ellps=WGS84), a quick plot shows this to not be the case.
<http://r-sig-geo.2731867.n2.nabble.com/file/n7587261/image822.png> 

Rather than being geographic coordinates, the coordinate values index to the
netCDF "lon" and "lat" values.
I can turn to ncdf4 to read the netCDF so as to read in the longitude and
latitude values for each indexed cell.  
The following code does this and casts all the data along with the actual
coordinates into a data frame, but it runs very slow.
Could there be a better way within the context of the raster package?

############# ncdf4 code that works, but is very slow (~16 hours of
processing) ##########
nc<-nc_open(filename=x)
nc
## has dimension time, j (384) cell index along second dimension, and i
(320) cell index along first dimension
dataLon <- ncvar_get( nc, "lon" )	# get the longitude lookup values
dim(dataLon)
dataLat <- ncvar_get( nc, "lat" )	# get the latitude lookup values
dim(dataLat)

###  Using nested loops, walk through the data frame inserting the lat and
lon values from the netCDF file
for(x in min(SICdf$xIndex):max(SICdf$xIndex)){  ## step through the time
steps and fill in the data
	for(y in min(SICdf$yIndex):max(SICdf$yIndex)){
		cat(x, y, fill=T)
		SICdf[SICdf$xIndex==x & SICdf$yIndex==y,]$lat<-dataLat[x,y][1]
		SICdf[SICdf$xIndex==x & SICdf$yIndex==y,]$lon<-dataLon[x,y][1]
	}
}
################################################################



-----
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650

AFischbach at usgs.gov
http://alaska.usgs.gov/science/biology/walrus
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/raster-working-with-netCDF-files-that-have-indexed-geography-tp7587261.html
Sent from the R-sig-geo mailing list archive at Nabble.com.