5 Operaciones

Tomado de Lovelace, Nowosad, and Muenchow (2019), capítulo 5.

Librerías usadas:

  • sp

5.1 Introducción

Usaremos la palabra geometría para denotar las características puramente espaciales, lo que significa que se ignoran los atributos (cualidades, propiedades de algo en una ubicación particular). Con ubicación/locación denotamos un punto, línea, polígono o celda de cuadrícula.

5.2 Unión (over)

La unión espacial ubica los índices o atributos de un objeto espacial en las ubicaciones de otro. En particular, el método over

over(x,y)
x %over% y

devuelve:

  1. En caso de que y no tenga atributos: los índices de y correspondientes a cada característica de x, o NA en caso de que no haya correspondencia;
  2. En caso de que y tenga atributos: un data.frame con los atributos de y correspondientes a las ubicaciones de x, o un registro NA en caso de que no haya correspondencia.

Correspondencia significa que dos características se cruzan espacialmente (tocar, superponer, cubrir, incluir, etc.).

over hace la siguiente operacion gráficamente en conjuntos:

drawing

Un uso típico de over es para seleccionar características de un objeto espacial que están sobre o dentro de otro.

Supongamos que queremos seleccionar todos los puntos en meuse que caen en meuse.grid, podríamos hacer esto de la siguiente manera:

library(sp)

r1 <- cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 
332618, 332413, 332349))
r2 <- cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683, 
331133, 331623, 332152, 332357, 332373))
r3 <- cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 
329783, 329665, 329720, 329933, 330478, 331062, 331086))
r4 <- cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791))

sr1 <- Polygons(list(Polygon(r1)),"r1")
sr2 <- Polygons(list(Polygon(r2)),"r2")
sr3 <- Polygons(list(Polygon(r3)),"r3")
sr4 <- Polygons(list(Polygon(r4)),"r4")
sr <- SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf <- SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2), 
    row.names=c("r1","r2","r3","r4")))

rm(meuse)
data(meuse)
coordinates(meuse) = ~x+y

sp::plot(meuse)
polygon(r1)
polygon(r2)
polygon(r3)
polygon(r4)

devuelve las concentraciones medias de metales pesados por polígono:

sp::over(sr, meuse[,1:4], fn = "mean")
##     cadmium   copper     lead     zinc
## r1 4.036000 44.48000 147.2600 475.8800
## r2 2.910526 40.22807 145.4035 452.0702
## r3 2.820833 36.08333 169.1667 484.2500
## r4       NA       NA       NA       NA

devuelve el número de puntos en cada polígono:

sapply(sp::over(sr, geometry(meuse), returnList = TRUE), length)
## r1 r2 r3 r4 
## 50 57 48  0
sp::over(sr, meuse)
##    cadmium copper lead zinc  elev       dist   om ffreq soil lime landuse dist.m
## r1    11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1      Ah     50
## r2     1.3     21   62  258 9.280 0.32057200  2.0     1    2    0      Ah    360
## r3     3.9     47  268  703 5.760 0.07568690  7.0     1    1    1      Ab     80
## r4      NA     NA   NA   NA    NA         NA   NA  <NA> <NA> <NA>    <NA>     NA

Se hizo: seleccionar los puntos de meuse que coinciden con (caen en) las celdas de sr.

También se puede trabajar con un SpatialGrid

rm(meuse.grid)
data("meuse.grid")
coordinates(meuse.grid) <- c("x","y")
sp::plot(meuse.grid)
sp::plot(sr, add = TRUE)

sp::over(sr, meuse.grid)
##    part.a part.b       dist soil ffreq
## r1      1      0 0.00135803    1     1
## r2      1      0 0.00000000    1     1
## r3      0      1 0.00241427    1     1
## r4     NA     NA         NA <NA>  <NA>
sp::over(sr, meuse.grid, fn = "mean")
##       part.a    part.b      dist soil ffreq
## r1 1.0000000 0.0000000 0.2233190   NA    NA
## r2 0.3853306 0.6146694 0.3881792   NA    NA
## r3 0.0000000 1.0000000 0.2224238   NA    NA
## r4        NA        NA        NA   NA    NA

5.3 Agregación espacial

Tomado de Pebesma (2021)

En el siguiente ejemplo, los valores de un grid con celdas de \(40m\times 40m\) se agregan a un grid con celdas de \(400m\times 400m\).

rm(meuse.grid)
data(meuse.grid)
#Convertimos a meuse.grid en SpatialPixelsDataFrame
gridded(meuse.grid) = ~x+y 

# Generamos la cuadricula para agregacion:
off <- gridparameters(meuse.grid)$cellcentre.offset + 20
gt <- GridTopology(off, c(400,400), c(8,11))
SG <- SpatialGrid(gt)

# Agregamos:
agg <- aggregate(meuse.grid[3], SG, mean)

Se muestran las gráficas del grid y los datos que se desea agregar:

sp::plot(meuse.grid[3])
sp::plot(SG,add = TRUE)

El agregado de funciones agrega su primer argumento sobre las geometrías del segundo argumento y devuelve una geometría con atributos. La función de agregación predeterminada (mean) se puede anular.

La siguiente muestra el resultado de esta agregación (agg, en colores) y los puntos (+) de la cuadrícula original (meuse.grid).

sp::plot(agg)
points(meuse,pch = 3)

5.4 Ejercicios

Importar los datos de Ambato urbano y las emergencias de tránsito reportadas en diciembre 2014. Con esta información:

library(rgdal) #readOGR
library(RColorBrewer)
library(classInt)
library(raster)


### Leo el shape:
uu <- "https://github.com/vmoprojs/DataLectures/raw/master/SpatialData/shape_Ambato_urbano.zip"
poligonos <- read_git_shp(uu)

uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/EmergenciasTransito.txt"
datos <- read.csv(uu, sep ="")
datos <- subset(datos,Month == 12)

  1. Contar el número de emergencias por cada parroquia urbana, debe obtener (verifique que el total sume 281):
##  ATOCHA-FICOA  LA PENINSULA     PISHILATA HUACHI LORETO CELIANO MONJE     LA MERCED 
##            21            15            37            37            60            27 
## SAN FRANCISCO     LA MATRIZ  HUACHI CHICO         Total 
##            19            31            34           281

  1. Generar una variable agregada por polígonos como una columna del data.frame asociado al polígono espacial.

  2. Transformar los polígonos a un objeto raster y generar un mapa de la variable Emergencias. Debe resultar el siguiente mapa: