[R] graphing netCDF files

[Paul Hiemstra]

Hi Steve,

If you read your netCDF files into R you end up with sp-classes which 
can be displayed using spplot. But you do not seem to use rgdal.

If you can make a data.frame with the x, y and z coordinates this can 
quite easily be transformed into an sp-class:

library(sp)
dat = data.frame(x = UTMx, y = UTMy, z = wat.data2001q1,,i])
coordinates(dat) = ~x+y   # "tell" spplot what the names of the columns 
with the x and y coordinates are
gridded(dat) = TRUE       # make clear it is a grid
spplot(dat)

For more details see the documentation for the sp-package, especially 
spplot. These kinds of questions are more suitable for the r-sig-geo 
mailing list and not the general r-help list.

hope this helps,

Paul

Steve_Friedman at nps.gov schreef:
> Hello
>
> I'm working with a large hydrological data set stored in a netCDF format.
> The file stores x and y coordinates in the UTM projected coordinate system,
> yet when I use image to graphically display the z variable, the image is
> distorted in the sense that it does not plot the map in the correct spatial
> organization.
>
> I'm wondering if I need to define the projection of the netCDF file with
> rgdal or proj4 routines first before I send it to the graphics device.
>   
Defining the projection is not needed
> My code is as follows:
>
>  q1_2001 <- open.ncdf("H:\\SKF_DESKTOP FILES\\My
> Documents\\EDEN\\EDEN\\Surfaces\\2000_q1.nc", readunlimi=FALSE) #opens ncdf
> file for reading
>    wat.data2000q1 <- get.var.ncdf(q1_2001,  verbose=FALSE ) # gets the real
> information
>
>  # GENERAL EXAMINATION OF HEADER DATA in the wat.data file
>    day <- get.var.ncdf(q1_2001, "time")   # length(day) 91 days in quarter
>    UTMx <-   get.var.ncdf(q1_2001, "x")   # columns (eastings)  # should
> return 405
>    UTMy <-   get.var.ncdf(q1_2001, "y")   # rows (northings)       #
> should return 287
>
> # plot first 91 days (3 months of the year)
>     for(i in 1:91) {
>        !is.na( image(UTMx, UTMy, z = wat.data2001q1[,,i], col=brewer.pal(8,
> "YlGnBu"),
>          axes=T, pty="s", ylab="UTM Northing", xlab="UTM Easting",
>          main = "First Quater 2001")  )
>          }
>
> As I indicated above the map is displayed on the graphics device. However
> the orientation is distorted pulling the x axis to wide and the y axis too
> tall.  How can I set the graphics device to know the orientation and
> scaling (if these are the correct terms) in order to display this map
> correctly?
>
> All insights will be greatly appreciated.
>
> Thanks
> Steve
>
> Steve Friedman Ph. D.
> Spatial Statistical Analyst
> Everglades and Dry Tortugas National Park
> 950 N Krome Ave (3rd Floor)
> Homestead, Florida 33034
>
> Office (305) 224 - 4282
> Steve_Friedman at nps.gov
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>   


-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:     +31302535773
Fax:    +31302531145
http://intamap.geo.uu.nl/~paul