mape_r-01: bivariate point pattern
Andree Valle-Campos
13/2/2020
1 objetivo
2 load packages
3 import data
3.1 explore the object
En tu computadora local, explora cuál es la salida de cada una de estas acciones
# always!
class(preston_crime)
# ppp is a class from the spatstat package
# explore with spatstat specific print
preston_crime %>% print.ppp()
# Get some summary information on the dataset
summary(preston_crime)
# list of attributes
attributes(preston_crime)
# explore attributes with $
preston_crime$markformat
# window is an exclusive element from ppp objects
preston_crime$window
# class it is?
class(preston_crime$window)
# object structure
str(preston_crime)4 transform data
# from ppp to tibble -----------------------
gpsrino <- preston_crime %>% as_tibble()
gpsrino %>% class()## [1] "tbl_df" "tbl" "data.frame"## # A tibble: 2,036 x 3
## x y marks
## <dbl> <dbl> <fct>
## 1 354558. 429074. Non-violent crime
## 2 353666. 429243. Non-violent crime
## 3 355801. 430015. Non-violent crime
## 4 352859. 429319. Non-violent crime
## 5 353560. 430174. Non-violent crime
## 6 351079. 431580. Non-violent crime
## 7 351046. 431497. Non-violent crime
## 8 354018. 429239. Non-violent crime
## 9 350119. 430782. Non-violent crime
## 10 355248. 431478. Non-violent crime
## # ... with 2,026 more rows# fromm tibble to sf -----------------------
library(sf)
house <- st_as_sf(gpsrino, coords = c("x", "y"), remove = F,
crs = 27561, agr = "constant")
# notes
# this always works for peru!
# peru: crs = 4326
# CRS("+proj=longlat +datum=WGS84")
house %>% class()## [1] "sf" "tbl_df" "tbl" "data.frame"## Simple feature collection with 2036 features and 3 fields
## Attribute-geometry relationship: 3 constant, 0 aggregate, 0 identity
## geometry type: POINT
## dimension: XY
## bbox: xmin: 350000 ymin: 430000 xmax: 360000 ymax: 430000
## epsg (SRID): 27561
## proj4string: +proj=lcc +lat_1=49.50000000000001 +lat_0=49.50000000000001 +lon_0=0 +k_0=0.999877341 +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs
## # A tibble: 2,036 x 4
## x y marks geometry
## <dbl> <dbl> <fct> <POINT [m]>
## 1 354558. 429074. Non-violent crime (354558 429074)
## 2 353666. 429243. Non-violent crime (353666 429243)
## 3 355801. 430015. Non-violent crime (355801 430015)
## 4 352859. 429319. Non-violent crime (352859 429319)
## 5 353560. 430174. Non-violent crime (353560 430174)
## 6 351079. 431580. Non-violent crime (351079 431580)
## 7 351046. 431497. Non-violent crime (351046 431497)
## 8 354018. 429239. Non-violent crime (354018 429239)
## 9 350119. 430782. Non-violent crime (350119 430782)
## 10 355248. 431478. Non-violent crime (355248 431478)
## # ... with 2,026 more rows5 exploratory plots
5.1 using geom_point
house %>%
#dplyr
mutate(marks=fct_relevel(marks,"Violent crime")) %>%
#ggplot2
ggplot(aes(x = x, y = y, color = marks, size = marks)) +
#geometry
geom_point() +
coord_fixed(ratio = 1) +
#aestetics
scale_color_manual(values = c("red","black")) +
scale_size_manual(values = c(1.5,0.5))
6 kernel smoothing
house %>%
#dplyr
mutate(marks=fct_relevel(marks,"Violent crime")) %>%
#ggplot2
ggplot() +
#geometry
geom_sf() +
coord_sf() +
facet_grid(~marks)
6.1 binwidth selection
6.1.1 histogram
binwidth =
The width of the bins. Can be specified as a numeric value or as a function that calculates width from unscaled x. Here, “unscaled x” refers to the original x values in the data, before application of any scale transformation. When specifying a function along with a grouping structure, the function will be called once per group. The default is to use the number of bins in bins, covering the range of the data. You should always override this value, exploring multiple widths to find the best to illustrate the stories in your data.
6.1.2 density
stat_density
Computes and draws kernel density estimate, which is a smoothed version of the histogram. This is a useful alternative to the histogram for continuous data that comes from an underlying smooth distribution.
bw =
The smoothing bandwidth to be used. If numeric, the standard deviation of the smoothing kernel. If character, a rule to choose the bandwidth, as listed in stats::bw.nrd().
# default: bw = "nrd0" method
# rule-of-thumb for choosing the bandwidth
# of a Gaussian kernel density estimator
house %>%
ggplot(aes(x = x)) +
stat_density(bw = 1000)
6.2 density 2d
stat_density_2d
Perform a 2D kernel density estimation using MASS::kde2d() and display the results with contours. This can be useful for dealing with overplotting. This is a 2d version of geom_density().
h =
Bandwidth (vector of length two). If NULL, estimated using MASS::bandwidth.nrd()
house %>%
ggplot() +
stat_density_2d(aes(x = x, y = y,fill = ..level..),
alpha=0.3,
geom = "polygon"
) +
coord_fixed(ratio = 1)
house %>%
ggplot() +
stat_density_2d(aes(x = x, y = y,fill = ..level..),
alpha=0.3,
geom = "polygon",
h = c(50,50)
) +
coord_fixed(ratio = 1)
house %>%
ggplot() +
stat_density_2d(aes(x = x, y = y,fill = ..level..),
alpha=0.3,
geom = "polygon",
h = c(5000,5000)
) +
coord_fixed(ratio = 1)
6.2.1 how to choose the better bandwidth?
# r function ----------------
ppp2sf <- function (ppp) {
data <- tibble::tibble(x = ppp$x, y = ppp$y)
if (!is.null(ppp$marks)) data$marks = ppp$marks
# data_sf <- st_as_sf(data, coords = c("x", "y"))
data_sf <- data
if (!is.null(ppp$window$bdry)) {
bnd <- as.matrix(as.data.frame(ppp$window$bdry[[1]]))
bnd <- rbind(bnd, bnd[1, ])
bnd <- st_sf(id = 1, geometry = st_sfc(st_polygon(list(bnd))))
} else {
bnd <- cbind(c(ppp$window$xrange, rev(ppp$window$xrange)),
rep(ppp$window$yrange, each = 2))
bnd <- rbind(bnd, bnd[1, ])
bnd <- st_sf(id = 1, geometry = st_sfc(st_polygon(list(bnd))))
}
return(list(data = data_sf, bnd = bnd))
}bw.scott
Use Scott’s rule of thumb to determine the smoothing bandwidth for the kernel estimation of point process intensity.
This function selects a bandwidth sigma for the kernel estimator of point process intensity computed by density.ppp.
# extract window to sf -----------
window_boundary <- ppp2sf(preston_crime) %>%
pluck(2) %>%
st_as_sf(remove = F,
crs = 27561, agr = "constant")
# determine the binwidth
library(maptools)
house_g <- house %>% select(geometry)
house_poly <- window_boundary %>% st_buffer(dist = 0) %>% st_union() #needs to be cleaner!
p.sp <- as(house_g, "Spatial") # Create Spatial* object
p.ppp <- as(p.sp, "ppp") # Create ppp object
Window(p.ppp) <- as.owin(as(house_poly, "Spatial"))
h_ppp <- bw.scott(p.ppp) #bw.ppl(p.ppp)
h_ppp## sigma.x sigma.y
## 472 422house %>%
ggplot(aes(x = x, y = y)) +
stat_density_2d(aes(fill = ..level..),
alpha=0.3,
geom = "polygon",
h = h_ppp
) +
coord_fixed(ratio = 1)
6.2.2 use multiple geometries
house %>%
ggplot(aes(x = x, y = y)) +
stat_density_2d(aes(fill = ..level..),
alpha=0.3,
geom = "polygon",
h = h_ppp
) +
geom_point() +
coord_fixed(ratio = 1)
house %>%
ggplot() +
stat_density_2d(aes(x = x, y = y,fill = ..level..),
alpha=0.3,
geom = "polygon",
h = h_ppp
) +
geom_sf(alpha=0.05,size=0.5) +
coord_sf()
house %>%
ggplot() +
stat_density_2d(aes(x = x, y = y,fill = ..level..),
alpha=0.3,
geom = "polygon",
h = h_ppp
) +
geom_sf(alpha=0.05,size=0.5) +
coord_sf() +
facet_grid(~marks)
7 bivariate point pattern
sigma
Standard deviation of isotropic smoothing kernel. Either a numerical value, or a function that computes an appropriate value of sigma.
# cargar paquete
library(readr)
library(spatstat)
# importar datos
preston_crime <- read_rds("pcrime-spatstat.rds.gz.rds")
# dividir ppp
crime_splits <- split(preston_crime)
crime_splits
plot(crime_splits)
# estimar densidad
crime_densities <- density(crime_splits,sigma=h_ppp)
crime_densities
plot(crime_densities)
# generar fracción de densidad ~ riesgo relativo
frac_violent_crime_density <- crime_densities[[2]] /
(crime_densities[[1]] + crime_densities[[2]])
#png("figure/fractional_density.png"),width = 600,height = 600)
plot(frac_violent_crime_density)
#dev.off()