[R] about netCDF

Agustin Lobo alobo at ija.csic.es
Fri Feb 22 10:54:40 CET 2002


After read.netCDF, use str() on the resulting
R object to see how it is organized. For example,
in one of my cases I have:

> str(npp2)
List of 9
 $ longitude   : num [, 1:90] -178 -174 -170 -166 -162 -158 -154 -150 -146
-142 ...
  ..- attr(*, "long_name")= chr "longitude"
  ..- attr(*, "units")= chr "degrees_east"
 $ latitude    : num [, 1:45] 88 84 80 76 72 68 64 60 56 52 ...
  ..- attr(*, "long_name")= chr "latitude"
  ..- attr(*, "units")= chr "degrees_north"
...
(output deleted here)
...
 $ npp         : num [1:5, 1:12, 1:45, 1:90] 9e+20 9e+20 9e+20 9e+20 9e+20
...
  ..- attr(*, "long_name")= chr "npp of carbon for each pft"
  ..- attr(*, "units")= chr "kg m-2 year-1"
  ..- attr(*, "missing_value")= num 9e+20
 $ npptot      : num [1:5, 1:45, 1:90] 9e+20 9e+20 9e+20 9e+20 9e+20 ...
  ..- attr(*, "long_name")= chr "total npp"
  ..- attr(*, "units")= chr "kg m-2 year-1"
  ..- attr(*, "missing_value")= num 9e+20
...
(more output deleted here)


There you can see, i.e., that npp is a numerical 4D array
with dimensions [1:5, 1:12, 1:45, 1:90], and that npptot 
is a 3D array with dimensions [1:5, 1:45, 1:90]

Then you can inquire more things on a
particular component, i.e.,
> attributes(npp2$npptot)
$dim
[1]  5 45 90

$"long_name"
[1] "total npp"

$units
[1] "kg m-2 year-1"

$"missing_value"
[1] 9e+20

Use R NA for missing values:

> range(a)
[1] -6.052019e-01  9.000000e+20
>  a[a == attributes(a)$"missing_value"] <- NA
> range(a,na.rm=T)
[1] -0.6052019  2.0153317

Add dimnames according to str(npp)

> dimnames(a) <- list(1901:1905,npp2$latitude,np2$longitude)

Global average Land npp for the period 1901-1905:
> mean(a,na.rm=T)
[1] 0.4340909

Average map of npp for  the period 1901-1905
a.medio <- apply(a,c(2,3),mean,na.rm=T)

1902 map of npp anomalies vs. the 1901-1905 mean:
a["1902",,] - a.medio

Plot a time profile:
plot(a[,"36","-6"],type="l",xlab="year",ylab="total npp (kg m-2 year-1)")

Plot a Latitude profile for year 1902 @ 22W:
plot(a["1902",,"22"],type="l",xlab="Latitude",ylab="Total Land NPP in 1902
Lon=22E (kg m-2 year-1)")

Plot a map of npp for 1903:
imagen(a["1903",,])

Check npp values interactively:
imaexplore(a["1903",,])

Hope this helps as an orientation. I can send you
npp2 if you need it.

Agus

Dr. Agustin Lobo
Instituto de Ciencias de la Tierra (CSIC)
Lluis Sole Sabaris s/n
08028 Barcelona SPAIN
tel 34 93409 5410
fax 34 93411 0012
alobo at ija.csic.es


On Thu, 21 Feb 2002, antonio wrote:

> Hi,
> 
> I would like to ask a couple of questions about netCDF package:
> 
> 1) I have COADS data in .cdf format. Data are from a 1ºx1º grid 
> for lat: x1-x2, lon: y1-y2, monthly values since 1960
> 2) I manage to open and read the file in R with your package without any 
> problem.
> 3) After opening and reading the file, how do I can manage the data. Is to 
> say, how do I can plot contours for an specific month or how do I can 
> average, for example, all the Jan, Feb, etc and then calculate anomalies. 
> 
> Thanks in advance,
> 
> Antonio Rodríguez
> CICEM Agua del Pino
> Huelva, Spain
> 
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> 

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list