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 [EMAIL PROTECTED] > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@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: [EMAIL PROTECTED] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo