[R-sig-Geo] (no subject)

[Edzer J. Pebesma]

Yong Li wrote:
> 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])
>   
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
> #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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>