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(inputim...@bbox[,1],
inputim...@grid@cellsize * SampleFactor,
inputim...@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)
outspatialgrids...@data <- samplingresultss...@data
outspatialgrids...@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(inputim...@data)
print(Incoming)
Outgoing = describe(samplingresultss...@data)
print(Outgoing)
browser()
#
dev.set(1) # plot in new window
grob <- spplot(OutSpatialGridSGDF, "band1",
col.regions=greys,cuts=length(greys))
plot(grob)
}
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
