[R-sig-Geo] (no subject)

Yong Li yong.li at unimelb.edu.au
Thu Apr 12 14:46:13 CEST 2007


require(maptools)

#read in shape file
pts   <- readShapePoints(system.file("shapes/baltim.shp",
package="maptools")[1])

#define the dimension of grids
x.min <- summary(pts)$bbox[1]
y.min <- summary(pts)$bbox[2]
x.max <- summary(pts)$bbox[3]
y.max <- summary(pts)$bbox[4]
cellsize <- 20
x.n   <- round((x.max-x.min)/cellsize)+1
y.n   <- round((y.max-y.min)/cellsize)+1

#produce polygons from grids
grd   <- GridTopology(c(x.min,y.min), c(cellsize,cellsize), c(x.n,y.n))
polys <- as.SpatialPolygons.GridTopology(grd)
plot(polys)
plot(pts, col="blue", add=TRUE )

index <- overlay(pts, polys[11])

#Then a. get the selected points in No.11 polygon,
      b. randomly just choose one point from above selected points,
      c. do the same for other polygons (rectangles), and
      d. put the finally selected points in another shape file.

The output would be that there is only one point in each rectangle. The
purpose of
doing this is to resample the points at different scales for
geo-statistics analysis.

Cheers

Yong




More information about the R-sig-Geo mailing list