[R-sig-Geo] Grids clipped with polygon

[Roger Bivand]

On Wed, 3 Aug 2005, Susumu Tanimura wrote:

> Hi there,
> 
> Grids clipped with polygon is just required for my kriging analysis,
> but I could not succeed so far.  I probed into gpclib, inside.owin() in
> spatstat, inpip() in splancs, symbolsInPoly() in maptools.
> 
> > library(maptools)
> > V.poly <- Map2poly(read.shape("vietnam-utm48n.shp"))
> Shapefile Type: Polygon   # of Shapes: 61
> > np <- sapply(V.poly, function(x) attr(x, "nPart"))
> > np
>  [1]   1 510   3   2   1   1   2   7   4   1   2   1   1   1   1   1   1 105   2
> [20]   5   1 682   1   1   1   2   7   8   3   2   3  46   1   3   8   7 138   1
> [39]   1   3   1   1   1   1   1   3   1   3  15   1   1   1   3   1   1  41   1
> [58]   2   1   1   1
> > try1 <- symbolsInPolys(V.poly, 5000)

Although you could use symbolsInPolys() for this, the function was written 
for display purposes. If you need to have a regular grid to cover the area 
of a shapefile of polygons, clipped to the polygon rings, I suggest using 
functions from the sp package:

> library(maptools)
Loading required package: foreign
> library(sp)
> x <- read.shape(system.file("shapes/sids.shp", package="maptools")[1])
Shapefile Type: Polygon   # of Shapes: 100 
> xSR <- as.SpatialRings.Shapes(x$Shapes, 
+ IDs=as.character(x$att.data$FIPS), proj4string=CRS(as.character(NA)))
> bbox(xSR)
         min       max
r1 -84.32385 -75.45698
r2  33.88199  36.58965
> grid <- GridTopology(cellcentre.offset=c(-85.4, 33.8), 
+ cellsize=c(0.1, 0.1), cells.dim=c(100, 40))
> grid.SP <- SpatialPoints(coordinates(grid), 
+ proj4string = CRS(as.character(NA)))
> summary(grid.SP)
Object of class SpatialPoints
Coordinates:
     min   max
s1 -85.4 -75.5
s2  33.8  37.7
Is projected: NA 
proj4string : [NA]
Number of points: 4000
> res <- overlay(grid.SP, xSR)
> summary(res)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   1.00   25.00   51.00   50.54   76.00  100.00 2734.00 
> plot(xSR)
> points(coordinates(grid.SP), pch=".", col="grey")
> points(coordinates(grid.SP)[!is.na(res),], pch=".", col="black", cex=2)

with coordinates(grid.SP)[!is.na(res),] being the matrix of grid 
coordinates within the polygons. Because GridTopology() lets you set the 
grid resolution directly, I think it will be better for non-display 
purposes than the function in the maptools package. By the way, gstat is 
using the classes and methods from the sp package in some cases, but this 
solution gives coordinates you can use generally. The sp package will also 
let you make SpatialPixels out of your points, which can next be filled in 
with missing data, and written out as an ArcInfo ASCII grid, if that is 
any help.

Best wishes,

Roger

> 
> Spatial Point Pattern Analysis Code in S-Plus
> 
>  Version 2 - Spatial and Space-Time analysis
> Error in if (counts[[i]] > 0) { : missing value where TRUE/FALSE needed
> In addition: There were 32 warnings (use warnings() to see them)
> 
> It is very happy if someone give me some advice about the error or
> alternative approach.  I just want to have grid trimmed by polygon
> boundary like meuse.grid in gstat.
> 
> --
> Susumu Tanimura
> Dept. of Socio-environmental Medicine
> Inst. of Tropical Medicine, Nagasaki University
> e-mail stanimura-ngs at umin.ac.jp
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
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