[R-sig-Geo] SP Spatial List operation - error in R 2.9 not in R 2.8 (download example)

[rick reeves]

Hello List:

While writing an R script to display and subsample a raster image,
I encountered an error in the Spatial_.xxxx classes that occurs in
R ver 2.9, but not in R ver 2.8.

Note: I installed both versions within the last 5 days, and ran
the update.packages() command after each installation.

Please download the script (and a small sample image .tif file) at:

   http://www.nceas.ucsb.edu/scicomp/SpSamplingGridQuest.zip

(complete script replicated at bottom of this message)

Here is the problem:

I developed this script using R 2.8, and the results were as expected:
A subsampled image. But when I ran the script under R 2.9.0, the
subsampled image is 'empty'. Checking the data objects, the first
statement to exhibit a problem is:

SamplingPoints <- as(SamplingGrid,"SpatialPoints")

For statement: < str(SamplingPoints) >

R version 2.8.1 produces this result:

Formal class 'SpatialPoints' [package "sp"] with 3 slots
  ..@ coords     : num [1:24021, 1:2] -1759592 -1758712 -1757832 
-1756952 -1756072 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] -1760032 5929 -1621872 140569
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr NA

But R version 2.9 produces this (incorrect) result. Compare the @coords slot
produced by each version:

Formal class 'SpatialPoints' [package "sp"] with 3 slots
  ..@ coords     : num [1:2, 1:2] -1759592 -1622312 6369 140129
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] -1760032 5929 -1621872 140569
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr NA
Warning message:
In data.row.names(row.names, rowsi, i) :
  some row.names duplicated: 
2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,
43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,.........

Have I made a coding mistake? Has anyone encountered a similar problem? 
Is there a fix underway?

Thanks for any information -

Rick R


-- 
Rick Reeves
Scientific Programmer/Analyst and Data Manager
National Center for Ecological Analysis and Synthesis
UC Santa Barbara
www.nceas.ucsb.edu
805 892 2533

The script: 

#########################################################################
# SubsampleImageFiles.r
#  
# This function reads a raster image file, and creates an output
# file with coarser spatial resolution. Technique used: Subsampling.
# Create a sampling grid at a different (usually coarser) resolution
# than the input image, 'overlay' the sampling grid on the input
# image, extract a point from the image at each grid location,
# and then write the sampling grid to an output file in .img format.
# 
# This example demonstrates the R feaures that read, write and display
# raster datasets AND use of the R Spatial objects used to store
# and manipulate raster images.
# 
# Note: as of May 22, 2009, this routine does not work properly
# under R version 2.9: 
# Statement 'SamplingPoints <- as(SamplingGrid,"SpatialPoints")': In R version 2.9,
#            creates incorrect Spatial Points object, with only 2 rows/columns. 
# 
# Author: Rick Reeves
# Date created: 29-Apr-2009                 
# Date modified: 20-May-2009                                                            
# NCEAS
#
#########################################################################
#
SubsampleImageFilesQuest <- function()
{
   library(rgdal)
   library(maptools)  
   library(Hmisc)
#
   SampleFactor <- 2 # A factor of 2 creates a new image with 1/2 the number of pixels as the input image
#
# Step 1) Read input image into a SpatialGridDataFrame object, then display.
# 
   InputImage <- readGDAL("NevadaSiteDEMAlbersSub.tif")      
#
# Create a basic 256-level grey scale 'ramp' for the DEM image
#   
   greys <- grey(0:256 / 256)   
   dev.set(1) # plot in new window      
   grob <- spplot(InputImage, "band1", col.regions=greys,cuts=length(greys))
   plot(grob)    
#
# Create a sampling grid, using the input image's spatial parameters  
#   
   SamplingGridTopology <- GridTopology(InputImage at bbox[,1],
                                       InputImage at grid@cellsize * SampleFactor,
                                       InputImage at grid@cells.dim / SampleFactor)
#
# Create the necessary data structures
#                                       
   SamplingGrid <- SpatialGrid(SamplingGridTopology)  
   OutSpatialGridSGDF <- as(SamplingGrid,"SpatialGridDataFrame")
   SamplingPoints <- as(SamplingGrid,"SpatialPoints")   # In R ver 2.9, this line generates a
#                                                       # SamplingPoints object with 2 rows, 2 cols. 
#                                                       # R ver 2.81 produces a  SamplingPoints object
                                                        # with same number of rows/cols as SamplingGrid
#                                                        
# Subsample the input image, create 
# a SpatialGridDataFrame object.
# 
   SamplingResultsSPDF <- overlay(InputImage,SamplingPoints)    
   OutSpatialGridSGDF at data <- SamplingResultsSPDF at data 
   OutSpatialGridSGDF at proj4string <- CRS(proj4string(InputImage))   
   gridded(OutSpatialGridSGDF) <- TRUE 
#
# use Hmisc library describe() function to compare the incoming and resampled images
#
   message("Summary of incoming and outgoing images (hit key to continue):") 
   Incoming = describe(InputImage at data)
   print(Incoming)
   Outgoing = describe(SamplingResultsSPDF at data)  
   print(Outgoing)
   browser() 
#   
   dev.set(1) # plot in new window
   grob <- spplot(OutSpatialGridSGDF, "band1", col.regions=greys,cuts=length(greys))
   plot(grob)   
}