Uso del paquete spatstat

spatstat

Se describe el uso del paquete spatstat para el uso de datos espaciales

Author

Affiliation

Yandira Cuvero

Escuela Politécnica Nacional

Published

August 14, 2023

Code

path <- '/Users/yandiracuvero/Documents/Pagina/'

En R

El paquete spatstat es un paquete colaborativo de R diseñado para realizar análisis de datos espaciales, el cual fue escrito por Adrian Baddeley y Rolf Turner. Este paquete es uno de los que cuenta con más contribuciones en R. Está formado por más de 1000 funciones y cuenta con un manual de 750 páginas. Para más detalle dirigirse a www.spatstat.org

La versión actual del paquete permite principalmente:

Crear, manipular y dibujar procesos puntuales.
Realizar análisis exploratorio a esos datos.
Simular procesos puntales bajo un cierto modelo.
Realizar pruebas de hipótesis.

Empezar

Antes de instalar spatstat es recomendable instalar el paquete deldir.

Code

# install.packages("deldir")
library(deldir)

Para empezar instalamos, cargamos la librería y revisamos los últimos avances del paquete.

Code

# install.packages("spatstat")
library(spatstat)

Podemos solicitar más detalles del paquete mediante el comando

Code

help(spatstat)

Finalmente, podemos solicitar información sobre las últimas actualizaciones del paquete

Code

#latest.news()

Demo

Para tener una idea de lo que se puede hacer con spatstat, podemos solicitar un demo del mismo. En el se muestran diferentes tipos de procesos puntuales, y cómo podemos estudiar sus propiedades.

Code

# demo(spatstat)

Procesos puntuales

Podemos trabajar con bases de datos que ya se encuentran incluidas dentro del paquete spatstat, o podemos crear procesos puntuales desde cero en base a bases de datos que tengamos a nuestra disposición.

Bases de datos disponibles

Para utilizar cualquiera de las bases de datos disponibles en el paquete, requeriremos solamente llamar a la base de datos mediante el comando data.

Code

library(spatstat)
data(cells)
plot(cells)

Las bases de datos disponibles se muestran a continuación

Crear nuevos objetos

Un objeto ppp que representa un Proceso Puntual dentro del paquete spatstat debe. tener las siguientes componentes:

ppp: patrón puntual plano
owin: Ventana de observación
tess: Partición considerada para el proceso puntual

Ventana

Podemos considerar dos tipos de ventanas, una regular (cuadrado, rectángulo, polígono, círculo, etc.) o una irregular (bordes políticos).

Para un cuadrado basta con indicar la longitud del mismo
Para una ventana que tenga forma de rectángulo, vamos a indicar los vértices del mismo, primero correspondiente al eje \(x\) y luego correspondiente al eje \(y\).
Para el círculo se debe indicar el diámetro y el centro del mismo.
Para un polígono se indica las lista de vértices del mismo, y se unirá con una recta los puntos entre ellos.

Opción 1: Definición mediante vértices

Code

cuadrado <- square(3)
rectangulo <- owin(c(0,20), c(0,10))
circulo <- disc(radius = 7, centre = c(0,5))
poligono <- owin(poly = list(x = c(2,5,7,3,0), y = c(0,0,3,6,3)))
  
par(mfrow = c(2, 2))
plot(cuadrado)
plot(circulo)
plot(rectangulo)
plot(poligono)

Opción 2: Mediante un shapefile

Code

library(maptools)
b <- readShapeSpatial(paste0(path, '/posts/2023-07-20-sf/data/nxparroquias.shp')) |>  subset(DPA_DESPAR=="QUITO")

bp  <- as(b, "SpatialPolygons")
bpp <- as(bp, "owin")
plot(bpp)

Datos

Opción 1: Cargar datos desde un .csv

Code

data(letterR)
data <- read.csv(paste0(path, '/posts/2023-08-14/base.csv'))
P <- ppp(data$x, data$y, xrange = c(0,12), yrange = c(0,12))
plot(P)

Opción 2: Desde un archivo .shp (Fuente: Geoportal.quito ).

Usamos el borde de la ciudad definido en la sección anterior. Corregimos la proyección para que tanto los puntos y el borde tengan la misma proyección. Luego, construimos el ppp al indicar las coordenadas del proceso puntual y la ventana que le correspondiente, con la ventana adecuada.

Code

library(maptools)
library(sp)
library(sf)

a_sf <- st_read(
      paste0(path,'/posts/2023-08-14/files/Paradas_DMQ.shp')
      )

Reading layer `Paradas_DMQ' from data source 
  `/Users/yandiracuvero/Documents/Pagina/posts/2023-08-14/files/Paradas_DMQ.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 3826 features and 9 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -78.58597 ymin: -0.379311 xmax: -78.392 ymax: -0.005249085
Geodetic CRS:  WGS 84

Code

target_crs <- st_crs("+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs")

a_sf_reprojected <- st_transform(a_sf, target_crs)
#print(st_crs(a_sf_reprojected))
class(a_sf_reprojected)

[1] "sf"         "data.frame"

Code

app <- ppp( x = st_coordinates(a_sf_reprojected)[,1],
            y = st_coordinates(a_sf_reprojected)[,2],
            marks = as.data.frame(a_sf)$Número,
            window = bpp)

Opción 3: A partir de una simulación de un modelo particular

Code

CRS <- rpoispp(100, owin = letterR)

Tenemos como resultado

Code

par(mfrow = c(1, 3))
plot(P, main = "Proceso puntual marcado" )
plot(app, main="Paradas de bus Quito", pch = 4)
plot(CRS, main = " CRS ")

Marcas

Si la base de datos original tiene la información de marca, se indica el nombre de la columna correspondiente, sino se puede agregar dentro de la creación del objeto ppp.

Opción 1: A partir de la base de datos directamente

Code

data <- read.csv(paste0(path, '/posts/2023-08-14/base.csv'))
M <- ppp(data$x, data$y, xrange = c(0,12), yrange = c(0,12), 
         marks = data$marca)

Opción 2: Asignando los valores de las marcas a un ppp conocido

Code

marks(P) <- rep(c("A", "B"), 10)

Tenemos como resultado

Code

par(mfrow = c(1, 2))
plot(M, main = "Desde datos")
plot(P, main = "Mediante función")

Unidades

Se puede agregar las unidades de los procesos puntuales, en ese caso, se debe indicar el nombre de las mismas en singular y plural.

Code

unitname(P) <- c("centímetro", "centímetros")
P

Marked planar point pattern: 20 points
marks are of storage type  'character'
window: rectangle = [0, 12] x [0, 12] centímetros

Partición

Una partición es una división del espacio de áreas disjuntas.

Existen varios tipos de partición que podemos realizar:

Política. Por ejemplo, la división de provincias, cantones, etc.
Artificial. Por ejemplo, una grilla de rectángulos.

Opción 1: Partición regular.

Code

tess(xgrid = 1:10, ygrid = 1:5)

Tessellation
Tiles are equal rectangles, of dimension 1 x 1 units 
9 by 4 grid of tiles
window: rectangle = [1, 10] x [1, 5] units

Code

quadrats(bpp, nx = 3, ny = 5)

Tessellation
Tiles are irregular polygons
13 tiles (irregular windows)
window: polygonal boundary
enclosing rectangle: [765008.8, 784851.8] x [9959531, 9994659] units

Opción 2: Otra opción es realizar una partición en función del proceso puntual, es decir, una región por cada punto.

Code

plot(dirichlet(cells), main = "Partición de Dirichlet o Voronoi")
points(cells, col = "blue")

Opción 3: Grafo de procesos puntuales.

Como vemos, no se forma una partición, sino un grafo.

Code

plot(delaunay(cells), main = "Partición de Delaunay")
points(cells, col = "blue")

Manejo de PPP

Por facilidad. para no tener que llamar a la base de datos por su nombre, podemos usar el comando

Code

X <- cells

Para “ver” las características de la base de datos, llamaremos poner un resumen de la base

Code

summary(X)

Planar point pattern:  42 points
Average intensity 42 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Window: rectangle = [0, 1] x [0, 1] units
Window area = 1 square unit

Podemos también dibujarla utilizando el comando plot integrado a esta base

Code

plot(X)

Fase exploratoria

Podemos contar el número de eventos para una grilla regular

Code

Q <-  quadratcount(X, nx = 5, ny = 4)
Q

            x
y            [0,0.2) [0.2,0.4) [0.4,0.6) [0.6,0.8) [0.8,1]
  [0.75,1]         1         2         1         3       1
  [0.5,0.75)       2         2         3         2       3
  [0.25,0.5)       2         3         3         2       2
  [0,0.25)         1         3         2         3       1

Ó, para una mejor comprensión podemos dibujar la cantidad de eventos en este

Code

plot(X)
plot(Q, add = TRUE, CEX = 2)

Intensidad

Si se sospecha que el proceso \(X\) es no homogéneo, la función de intensidad o la medida de intensidad pueden ser estimados mediante técnicas no paramétricas.

Conteo de cuadrantes

En este caso, la ventana \(W\) se divide en subregiones \(B_1, \ldots, B_n\) iguales al área. Se conteo el número de eventos en cada cuadrante,

\[ n_j = n(x \cap B_j) \qquad \text{ para } j = 1, \ldots, m. \]

este es un estimador insesgado de los valores \(\Lambda(B_j)\).

Ejemplo

Code

Q <- quadratcount(cells, nx = 6, ny = 5)
plot(cells, cex = 1, pch = "+")
plot(Q, add = TRUE, cex = 1, col = "blue")

Kernel

La estimación usual del kernel de la función de intensidad es

\[ \tilde{\lambda} = e(u) \sum_{i=1}^n \kappa(u - x_i), \]

donde \(\kappa(u)\) es el kernel (una función de densidad arbitraria, y

\[ e^{-1}(u) = \int_W \kappa(v-u) dv \]

es el efecto del sesgo en el borde. Claramente, \(\tilde{\lambda}(u)\) es un estimador insesgado de

\[ \lambda^{*}(u) = e(u) \int_W \kappa(u-v) \lambda(v) dv \]

que es una versión suavizada de la verdadera función de intensidad \(\lambda(u)\). La elección de suavizar el kernel tiene como desventaja un incremento en el sesgo y/o en la varianza.

Ejemplo

Por default, se considera el kernel Guassiano isotrópico.

Code

den <- density(cells, sigma = 70)

par(mfrow = c(1, 3))
plot(den, main = "Gráfico de pixel")
plot(cells, add=TRUE, cex = 1, pch = "+")

contour(den, main="Gráfico de contorno")

persp(den, main="Gráfico de nivel")

Función de Ripley

Code

K <- Kest(X)
plot(K)

Otros paquetes a explorar

Referencia

Baddeley, A. (2008). Analysing spatial point patterns in R. Technical report, CSIRO, 2010. Version 4. Available at www.csiro.au/resources/pf16h.html