[R-sig-Geo] SpatialGridDataFrame Help

[Roger Bivand]

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.Bivand at nhh.no