Skip to content

Overlay of SpatialGrid/Raster and polygons

7 messages · Ned Horning, Robert J. Hijmans, Alan Swanson +2 more

Original
Ned Horning
# 
Greetings,

I am looking for a way to use R to extract pixel values (currently in a 
large RasterStack object) that fall under polygons (currently a 
SpatialPolygonsDataFrame object). I seem to recall a discussion about 
using overlay to do this but I can't find a method that would work.

Any insight would be appreciated.

Ned
Robert J. Hijmans
# 
Hi Ned,

Here is an approach to get values from a RasterStack to all cells in
each polygon:

library(raster)
# a polygon
data(meuse.riv)
pol <- SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)), "x")))
# a raster
r1 <- raster(system.file("external/test.ag", package="sp"))
r2 <- sqrt(r1)
# a stack
s <- stack(r1, r2)

# convert the polygon to a RasterLayer
rr <- polygonsToRaster(pol, s)

par(mfrow=c(1,2))
plot(r1)
plot(pol, add=TRUE)
plot(rr)

# extract points that are not NA
pts <- rasterToPoints(rr) # see additional arguments to select a
subset (useful for very large rasters)
# perhaps subsample your points
sampl <- sample(1:length(pts[,1]), min(100, length(pts[,1])))
pts <- pts[sampl,]

v <- xyValues(s, pts[,1:2])
# if you have mutiple polygons, bind the polygon ID to the raster values:
v <- cbind(pts[,3], v)

# You could also sample points with spsample
pts <- spsample(pol, 100,  "random")
# remove duplicate cells
cells <- unique(cellFromXY(s, pts))
v2 <- cellValues(s, cells)

# or sample from a SpGDF if you can create it from the RasterLayer (if
it is not too big)
spdf <- as(rr, 'SpatialGridDataFrame')
pts <- spsample(spdf, 100,  "random")
v3 <- xyValues(s, pts)


For extremely large rasters, rasterToPoints could fail. If so, first
create a RasterLayer with random values with calc and "fun=runif" and
make values F for e.g x > 0.01, overlay that with rr, and try
rasterToPoints again...

Hth, Robert
On Mon, Nov 23, 2009 at 12:16 PM, Ned Horning <horning at amnh.org> wrote:
Alan Swanson
# 
Ned, here is how I have done that using the overlay function. Not sure 
if it would work with a rasterstack. 
Alan

xy.locs <- readOGR(poly.shapefile.name,drop_unsupported_fields=T)
img<-readGDAL(image.name)
y<-slot(img,"data")
n.bands <- ncol(y)
bnames <- paste(varname,"_b",1:n.bands,sep="") }
out.frame <- 
as.data.frame(matrix(NA,nrow=nrow(xy.locs),ncol=n.bands,dimnames=list(NULL,bnames)))   

for(m in 1:(n.poly <- length(slot(xy.locs,"polygons")))){
    gc()
    z<-overlay(img,xy.locs[m,])
    for(j in 1:n.bands){
        xx <- mean(y[z==1 & !is.na(z),j],na.rm=T)
        out.frame[m,j] <- xx
        } # end inner loop over bands #
    } # end loop over individual polygons #
Ned Horning wrote:
Zia Ahmed
# 
I  am trying to install "raster" package, but  it is not available  in 
this site.  http://R-Forge.R-project.org .

Zia
Robert J. Hijmans wrote:
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20091123/56408ae6/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: zua3.vcf
Type: text/x-vcard
Size: 281 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20091123/56408ae6/attachment.vcf>
Robert J. Hijmans
# 
ZIa,

You can download the source code and install yourself here:
http://r-forge.r-project.org/R/?group_id=294

Or use this: install.packages("raster", repos="http://R-Forge.R-project.org")

But I think you need to have the latest version of R (currently 2.10)
installed for that to work

Robert
On Mon, Nov 23, 2009 at 2:08 PM, Zia Ahmed <zua3 at cornell.edu> wrote:
Zia Ahmed
# 
Hi all,

I want project several  several raster(.tif)  and vector (.shp)  files 
in R using rgdal driver. One set of  files are in  geographic coordinate 
system (WGS84) other are in LCC-Bangladesh coordinate system.

.I want to re-project them to EPSG:3104 (Gulshan 303/ BTM).  If some 
help me this regard it will appreciated.

Thanks

Zia

Below I mentioned parameter of LCC-Bangladesh coordinate system:
Projection: Lambert_Conformal_Conic
False_Easting: 2743186.000000
False_Northing: 914395.000000
Central_Meridian: 90.000000
Standard_Parallel_1: 23.176944
Standard_Parallel_2: 28.822770
Latitude_Of_Origin: 26.000000
Linear Unit: Meter (1.000000)

Geographic Coordinate System: GCS_Everest_Bangladesh
Angular Unit: Degree (0.017453292519943299)
Prime Meridian: Greenwich (0.000000000000000000)
Datum: D_Everest_Bangladesh
Spheroid: Everest_Adjustment_1937
 Semimajor Axis: 6377276.344999999700000000
  Semiminor Axis: 6356075.413140240100000000
  Inverse Flattening: 300.801699999999980000



-------------- next part --------------
A non-text attachment was scrubbed...
Name: zua3.vcf
Type: text/x-vcard
Size: 281 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20091123/b8a8f59e/attachment.vcf>
Roger Bivand
#
On Mon, 23 Nov 2009, Zia Ahmed wrote:

            

      
      
      
    

        
For the vector files, make sure that the Spatial*DataFrame objects read 
from the shapefiles have valid coordinate reference system descriptions. 
Then use spTransform() in rgdal to CRS("+init=epsg:3104"). Note that I see 
3106 not 3104 in EPSG. Further, EPSG does not define a datum, so more work 
will be required to establish the actual input and output CRS values.

See http://www.asprs.org/resources/GRIDS/ for March 2008, which includes 
values for a three-parameter +towgs84= transformation for Gulshan (watch 
the signs!). From Cliff Mugnier's description, it looks as though your LCC 
below is "+init=epsg=24375", # Kalianpur 1937 / India zone IIb, which does 
not include a +towgs84= conversion - more searching will be needed. There 
is a +towgs84= for Kalianpur 1975, but the ellipsoid is different from 
yours. I see TOWGS84[214.0, 804.0, 268.0, 0.0, 0.0, 0.0, 0.0] in 
http://apps.who.int/tools/geoserver/srsHelp.do, but you need someone with 
knowledge rather than Google!

Typically, you need to find a ground control point on your vector map, 
transform it to "+proj=longlat +datum=WGS84", and use writeOGR() with the 
"KML" driver to view in Google Earth. Because you can zoom in, you should 
be able to see when the coordinate reference system is adequate. See also 
http://spatialreference.org/.

You cannot project raster data, you have to warp it to a regular grid in 
the target coordinate reference system. This involves interpolation, so 
see for example projectRaster() in the raster package on R-Forge for 
ideas.