Tuesday, 10 October 2017

raster - how to extract a variable at specific time and depth from netcdf in r



I have the following netcdf file containing 7 variables and 4 dimensions.


[1] "file data.nc has 4 dimensions:"
[1] "depth Size: 12"
[1] "time Size: 37"
[1] "x Size: 305"
[1] "y Size: 432"
[1] "------------------------"
[1] "file data.nc has 7 variables:"
[1] "int stereographic[] Longname:stereographic Missval:NA"
[1] "float longitude[x,y] Longname:longitude Missval:1e+30"

[1] "float latitude[x,y] Longname:latitude Missval:1e+30"
[1] "short chla[x,y,depth,time] Longname:chla Missval:-32767"
[1] "short nitrat[x,y,depth,time] Longname:nitrat Missval:-32767"
[1] "short phosphat[x,y,depth,time] Longname:phosphat Missval:-32767"
[1] "short pbiomass[x,y,depth,time] Longname:pbiomass Missval:-32767"

I would like to extract several months at a specific depth and convert them to a ASCII raster. I am able to use this command myraster<-raster("data.nc", varname="nitrat", lvar=1, band=2, level=20) but I am confused about which depth and time I am actually extracting, it is unclear to me how to use the parameters of the raster function to point to the specific time and depth. Also I am not able to read a variable using the function get_att and I get this message


test<-open.ncdf("data.nc")
Error in ncatt_get(test, "time") :
Error, first passed argument must be an object of class ncdf4


I would like to be able to get a list of the different time and depth so I can use the correct parameters in the raster function.I have both "ncdf","ncdf4" packages installed.I am working in Windows with R 3.0



Answer



It might be easier to see what is going on when you create a RasterBrick (multiple layers), either representing depth, or time. The number of layers should make it clear which is is. the z-values too.


You could start with


b <- brick("data.nc", varname="nitrat")

and go from there.


The error in hour ncdf code is because you mix up the ncdf package (open.ncdf) with the ncdf4 package (ncatt_get). You need to pick one and stick with it.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...