On Wed, 25 Aug 2010, dkod wrote:
I am trying to develop simple tools to download and plot climate data based
on global grid of 2 degrees. Here's a link to the source image
http://data.giss.nasa.gov/work/gistemp/NMAPS/tmp_GHCN_GISS_HR2SST_1200km_Anom07_2010_2010_1951_1980/GHCN_GISS_HR2SST_1200km_Anom07_2010_2010_1951_1980.gif
link
NASA provides a 5 column text file for each image of monthly global
temperature anomalies. I have downloaded a sample data file that can be
viewed http://processtrends.com/Files/global_lota_map_07_36.txt here .
I have been using Bivand et al's Applied Spatial Data Analysis With R to try
to speed up my learning curve, however, I am stumped and need some help.
Here's my script so far:
#################################################################
### Read GISS global temp anom for month by 2 degree grid
library(fields);library(mapproj)
library(sp); library(rgdal)
# Step 1: Read source data file
link <- "http://processtrends.com/Files/global_lota_map_07_36.txt";
rdf <- read.table(link, skip = 1, header=T)
a_rdf <- data.frame(rdf[,c(1,2,5)])
names(a_rdf) <- c("i", "j","anom")
a_rdf$anom[a_rdf$anom==9999.0000] <- NA # convert all 9999.0000 to
NA
## Step 2: Setup sp() classes for GridTopology & SpatialGrid
grd <- GridTopology(cellcentre.offset=c(-179,-89), cellsize=c(2,2),
cells.dim=c(180,90))
gr_sp <- SpatialGrid(grid=grd, proj4string = CRS(as.character(NA)))
a_rdf_sp <- SpatialGridDataFrame(grd, a_rdf)
## Step 4: Create SpatialGridDataFrame
fullgrid(a_rdf_sp) <- T
slot(a_rdf_sp, "grid")
## Step 5: Generate Plot
image.plot(a_rdf_sp["anom"])
##############################################################
It runs without error mesages until the image.plot() command. here's the
error message I get
> image.plot(a_rdf_sp["anom"])
Error in if (del == 0 && to == 0) return(to) :
missing value where TRUE/FALSE needed
>
Well, either use the image() method in sp for this class of object, since
you have taken trouble to create it, or coerce the object to the input
required by image.plot() in the fields package:
library(sp)
data(meuse.grid)
coordinates(meuse.grid) <- c("x", "y")
gridded(meuse.grid) <- TRUE
fullgrid(meuse.grid) <- TRUE
mg1 <- as.image.SpatialGridDataFrame(meuse.grid["dist"])
library(fields)
image.plot(mg1)
This has relevance for another current thread on conversion between
classes. If package authors and maintainers do not wish to provide
coercion to or from sp, it would not be in the spirit of the R project to
jump in boots first. Documentation of possible routes is of essence here,
and could happily be done on the R Wiki site under spatial data - for
instance, an interested person could post the code above as an answer to
the sp -> fields question for this class and that function. Agreed, the
routes that exist could also be made a little more consistent, but since
not all the targets for conversion are S3 or S4 classes ("image" is not a
class, the output is a list with x, y, z components), it is hard to
find a consistent definition of consistent!
The relationship between spatstat and sp can serve as an example of how to
handle things pragmatically and with respect for the value of different
ways to represent data; coercion methods have been added as needed to
maptools, and work satisfactorily.
Data representations do differ, and there are often good, substantive
reasons for the differences, based on different literatures in different
disciplines or traditions. Giving users the impression that this is not
the case blinds them to the potential misconceptions involved (in
GIScience - ontological mismatch). So exposing users to the discomfort of
having to think through which representations may differ does have a real
motivation - say in drilling down to say matrix <-> graph representations,
as in Adrian Baddeley's post this morning, where thinking through
representations yielded fruitful insights.
Roger
I'm not clear how to best establish a SpatialGridDataFrame for global data
series from 2 or 5 degree files. I'd like to be able to plot using mercator
or other projection.
I'd also like to be able to add land forms.
I'd appreciate any help in correcting/improving this script. My goal is a
clear, easy to follow R script that can be used by R users to download and
work with global climate data files from NASA, NOAA and other agencies.
D Kelly O'Day
http://chartsgraphs.wordpress.com
http://processtrends.com
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: roger.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
