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

[Tyler Frazier]

Hi Tony: Previously I used layer1 <- readGDAL(“*myfile.nc") then <- raster() and then <- stack(layer(length(1:n)) but then it turned out stack() should have worked just fine to being with and I had my syntax a bit off.  I’m not sure in your case, but crs(my_file) <- “+proj…” could be of use.  Also, are you concatenating a bunch of files from ESGF after the wget?  Maybe something is going on with the original file.  16 hours seems like an awful long time to bring the file into R using ncdf4::  HTH , Tyler

On 10 Oct 2014, at 11:32 pm, Anthony Fischbach <afischbach at usgs.gov> wrote:

> 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.
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo