-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRpoints.Rmd
112 lines (91 loc) · 4.05 KB
/
Rpoints.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
---
title: "Generate Random Sampling Points by Minimum Distance"
author: "Anthony Caravaggi"
date: "`r format(Sys.Date())`"
output: github_document
---
Libraries
```{r, message = FALSE}
devtools::install_github("dkahle/ggmap") # v2.7 required
library(ggmap)
library(spatstat)
library(rgdal)
library(rgeos)
library(raster)
library(maptools)
library(assertthat)
```
This demonstration will go through the generation of random points separated by a minimum distance threshold, step-by-step. A function, Rpoints, condensing the process can be found at the bottom of the page.
Create bounding polygon.
```{r}
x <- c(-3.541238, -3.547341, -3.561397, -3.572132, -3.572495, -3.551687, -3.566528, -3.599789, -3.592162,
-3.604156, -3.629145, -3.618190, -3.585571, -3.572779, -3.525398, -3.504483, -3.525454, -3.541238)
y <- c(51.68106, 51.69050, 51.7003, 51.69617, 51.68800, 51.67405, 51.64949, 51.65832, 51.67560,
51.68177, 51.71373, 51.72463, 51.73346, 51.73788, 51.71745, 51.69461, 51.67548, 51.68106)
xy <- cbind(x,y)
p1 <- Polygon(xy)
poly <- SpatialPolygons(list(Polygons(list(p1), ID = "a")), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
```
```{r, echo = FALSE}
plot(poly)
```
Generate random points with a minimum distance between each point. The rSSI function uses an inhibition distance based on a Simple Sequential Inhibition point process. I used it here for simplicity; other limited point-generation methods could be used instead.
```{r}
samp1 <- rSSI(0.0025, 10, poly)
samp1 <- cbind(samp1$x, samp1$y) # Extract coordinates
samp1 <- SpatialPoints(samp1, crs(poly)) # Transform to SpatialPoints shapefile
```
```{r echo = FALSE}
plot(poly)
plot(samp1, add = TRUE)
```
Check that the minumim threshold has not been violated.
```{r}
spDists(samp1)*1000
```
```{r echo = FALSE}
spDists(samp1)*1000
```
Extract point coordinates and join to SpatialPoints object to create SpatialPointsDataFrame
```{r}
df <- data.frame(x = samp1@coords[,1], y = samp1@coords[,2])
samp1 <- SpatialPointsDataFrame(samp1, df)
```
We might as well go ahead and map these to a Stamen Map image. This demonstration previously mapped on top of a Google Maps image but Google stopped providing free API access in July 2018. If you want to use this method wit Google Maps then you'll require an API key and to use the `get_googlemap` function. Create a bounding box, download map tiles image based the bounding box at a given resolution (zoom), and plot polygon and points.
```{r eval = FALSE}
bbx <- c(left=-3.664351,bottom=51.631480,right=-3.467390,top=51.747602)
tre <- get_stamenmap(bbx, zoom =12, maptype = "watercolor")
sDat <- data.frame(samp1@coords)
ggmap(tre) + geom_polygon(aes(x=x, y=y), data=poly, fill="red", alpha=.5) +
geom_point(aes(x = sDat[,1], y = sDat[,2]), data = sDat, col="orange")
```
```{r echo = FALSE}
bbx <- c(left=-3.664351,bottom=51.631480,right=-3.467390,top=51.747602)
tre <- get_map(bbx, source = "stamen", zoom =12, maptype = "terrain")
sDat <- data.frame(samp1@coords)
ggmap(tre) + geom_polygon(aes(x=x, y=y), data=poly, fill="red", alpha=.5) +
geom_point(aes(x = sDat[,1], y = sDat[,2]), data = sDat, col="orange")
```
These steps - minus the mapping - have been condensed into one function, Rpoints, below. It includes an additional step - the extraction of the lowest non-zero vaue from the distance matrix, which is then added to the dataframe.
```{r eval}
# Function to generate random points within a polygon buffered by a minimum distance
# Attaches a dataframe of distances-between-points
#
# d = inherited distance (see rSSI documentation)
# n = number of points
# p = polygon
#
# E.g.
# spat <- Rpoints(d = 0.05, n = 50, p = poly1)
Rpoints <- function(d = 0.0025, n = 10, p){
s <- rSSI(d, n, p)
y <- SpatialPoints(cbind(s$x, s$y), crs(p))
t <- data.frame(spDists(y)*1000)
d <- data.frame(dist = apply(t, 1, function (v) min(v[v > 0])), # Collapse distances to non-zero minimum
x = y@coords[,1],
y = y@coords[,2])
o <- SpatialPointsDataFrame(y,d)
}
spat <- Rpoints(d = 0.001, n = 30, p = poly)
head(spat@data)
```