Skip to content

SpatialGridDataFrame Help

6 messages · Roger Bivand, dkod

Original
dkod
# 
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
  > 

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
#
On Wed, 25 Aug 2010, dkod wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      dkod
    

    

      
      #
    

  


  

      
Roger

Thank you for your quick and helpful response.

I now have image/plot() working. I need some help in how to go from
rectangular grid representation to a mercator (or other) projection. Since
my source grid is 2 degree long/lat based, I'd like to use a projection. 


How do I get a proper projection? Here's the line I used:

 gr_sp <- SpatialGrid(grid=grd, proj4string = CRS(as.character(NA))) 

Sorry for such a basic question.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 26 Aug 2010, dkod wrote:

        

            

      
      
      
    

        
No, much more work needed. If the input is in decimal degrees, 
geographical coordinates, and a known datum, you must warp (interpolate) 
to a planar projection. You might do this by projecting the input as a 
SpatialPointsDataFrame from +proj=longlat and known +datum to your target 
projection, then interpolate using your chosen interpolator to a grid 
defined in the target projection. A grid in geographical coordinates will 
effectively never be a grid in projected coordinates.

Roger

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      dkod
    

    

      
      #
    

  


  

      
Roger

Could you give me a hint on how I can add a landform background to my grid
plot? 

How do I use the maps packages'   world coordinates to show the landforms.

Thanks again

Kelly O'Day

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 26 Aug 2010, dkod wrote:

        

            

      
      
      
    

        
See map2SpatialLines() in maptools, work out whether to transform, or at 
least to assign the correct CRS. Note that image.plot() - if you are still 
using that function - does not stretch the Northings (y-axis), so if you 
are away from the Equator, your output will look flattened. Plotting 
functions in the sp and maps packages do stretch when the data are known 
to be in geographical coordinates.

Roger