Skip to content

(no subject)

2 messages · Yong Li, Edzer Pebesma

Original
Yong Li
# 
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
Edzer Pebesma
#
Yong Li wrote:
This should give you the indices of pts that fall in poly[11].
You can randomly sample this by
one.index = sample(index, 1)
and select the point by
pts[c(one.index,two.index),]
etc.

HTH,
--
Edzer