3 Objetos espaciales

Tomado de R. S. Bivand et al. (2013), capítulo 2

Librerías usadas:

  • sp
  • maps
  • maptools
  • raster

3.1 Introducción

Topógrafos, navegantes e ingenieros militares y civiles refinaron los conceptos fundamentales de la geografía matemática, establecidos a menudo hace siglos por algunos de los fundadores de la ciencia.

El uso de datos espaciales como vehículo comercial se ha visto impulsado a principios del siglo actual por la penetración de la banda ancha móvil (y cable) de los consumidores y las granjas de servidores distribuidos, con ejemplos como Google Earth ™, GoogleMaps ™ y otros.

Además, el espacio y otras tecnologías aerotransportadas han aumentado enormemente los volúmenes y tipos de datos espaciales disponibles:

Las fuentes de información entregan puntos espaciales en un sistema de referencia de coordenadas:

  • Combinando pares de puntos se forman segmentos de líneas.
  • Combinando segmentos de líneas se forman polí-líneas, redes, polígonos, cuadrículas.
  • Las cuadrículas se pueden definir dentro de un polígono regular, generalmente un rectángulo, con una resolución determinada: el tamaño de las celdas de la cuadrícula.

Sobre estas definiciones se conforman los modelos de datos en GIS donde diferentes elecciones implican compensaciones entre precisión, viabilidad y costo.

Cuando hay varias observaciones sobre el mismo atributo a diferentes alturas o profundidades, con frecuencia se tratan como capas separadas.

3.2 Clases y métodos en R

Pese a lucir como una calculadora, R tiene funciones y operadores detrás:

pi*10^2 # como lo escribimos
## [1] 314.1593
"*"(pi, "^"(10, 2)) # lo que sucede detrás
## [1] 314.1593

Una de las abstracciones fundamentales utilizadas en R es la formula.

Una fórmula suele tener dos caras, con una variable de respuesta a la izquierda del operador ~ (tilde) y, en este caso, una variable determinante a la derecha:

class(dist ~ speed)
## [1] "formula"

Se puede usar en contextos como modelización o:

plot(dist ~ speed, data = cars)

La ventaja central de las nuevas clases (S4) es que tienen definiciones formales que especifican el nombre y el tipo de los componentes, llamados slots, que contienen. El paquete sp usa clases S4 (en su mayoría) para representar objetos espaciales.

3.3 Objetos espaciales

La clase fundamental es Spatial con dos slots:

  • El primero es una matriz con nombres de columnas c("min", "max") con al menos dos filas, eje \(x\) y eje \(y\) respectivamente.

  • El segundo es un objeto de clase CRS que define el sistema de referencia de coordenadas y puede ser NA.

  • Operaciones con objetos de clase Spatial* deben actualizarse o copiar estos valores en e nuevo objeto Spatial*.

Se puede usar getClass para obtener la información del objeto espacial y sus subclases:

# library(sp)
getClass("Spatial")
## Class "Spatial" [package "sp"]
## 
## Slots:
##                               
## Name:         bbox proj4string
## Class:      matrix         CRS
## 
## Known Subclasses: 
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2

Veamos brevemente la clase CRS:

getClass("CRS")
## Class "CRS" [package "sp"]
## 
## Slots:
##                 
## Name:   projargs
## Class: character
  • La clase tiene una cadena de caracteres como único valor en el slot, que puede ser un valor NA. Para las coordenadas geográficas, la string más simple es + proj = longlat, con longlat. Veamos un ejemplo:
m <- matrix(c(0, 0, 1, 1), ncol = 2, dimnames = list(NULL, c("min", "max")))
crs <- sp::CRS(projargs = as.character(NA))
crs
## CRS arguments: NA
S <- sp::Spatial(bbox = m, proj4string = crs)
S
## class       : Spatial 
## features    : 1 
## extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs         : NA

Si se sabe que el objeto es proyectado:

bb <- matrix(c(320, 85, 380, 95), ncol = 2, dimnames = list(NULL, c("min", "max")))
sp::Spatial(bb, proj4string = sp::CRS("+proj=longlat"))
## class       : Spatial 
## features    : 1 
## extent      : 320, 380, 85, 95  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs

La definición de la clase Spatial no incluye simbología, tamaño, formas o colores. En sp se maneja estos temas cuando se despliegan los objetos gráfricamente.

3.4 Puntos espaciales

Un punto bidimensional puede describirse mediante un par de números (x, y), definidos sobre una región conocida:

un punto es un objeto geométrico de dimensión 0 y representa una única ubicación en el espacio de coordenadas (Herring and others (2011))

[VIDEO] LATITUD Y LONGITUD - ¿PORQUÉ ESTÁN, EN GRADOS?

  • Para representar fenómenos geográficos, la región máxima conocida es la superficie de la tierra, y el par de números medidos en grados son una coordenada geográfica, que muestra dónde está nuestro punto en el globo.

  • Las coordenadas geográficas pueden extenderse desde la latitud 90º a -90º n la dirección norte-sur y desde la longitud 0º a 360º o, de manera equivalente, desde -180º a 180º en la dirección este-oeste.

Usando la función estándar read.table, leemos en un archivo de datos con las posiciones de los espejos CRAN en todo el mundo en 2005. Extraemos las dos columnas con los valores de longitud y latitud en una matriz:

uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/CRAN051001a.txt"
CRAN_df <- read.table(url(uu), header = TRUE)
CRAN_mat <- cbind(CRAN_df$long, CRAN_df$lat)
row.names(CRAN_mat) <- 1:nrow(CRAN_mat)
str(CRAN_mat)
##  num [1:54, 1:2] 153 145 16.3 -49.3 -42.9 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:54] "1" "2" "3" "4" ...
##   ..$ : NULL

Recordemos que el objeto Spatial tenía dos slots, una baja que indica los límites y un CRS, SpatialPoints añade un slot llamado coords donde se puede insertar la matriz de coordenadas.

getClass("SpatialPoints")
## Class "SpatialPoints" [package "sp"]
## 
## Slots:
##                                           
## Name:       coords        bbox proj4string
## Class:      matrix      matrix         CRS
## 
## Extends: "Spatial", "SpatialVector", "SpatialPointsNULL"
## 
## Known Subclasses: 
## Class "SpatialPointsDataFrame", directly
## Class "SpatialPixels", directly
## Class "SpatialPixelsDataFrame", by class "SpatialPixels", distance 2

Un objeto SpatialPoints puede contener solo un punto.

llCRS <- sp::CRS("+proj=longlat +ellps=WGS84")
CRAN_sp <- sp::SpatialPoints(CRAN_mat, proj4string = llCRS)
sp::summary(CRAN_sp)
## Object of class SpatialPoints
## Coordinates:
##                  min      max
## coords.x1 -122.95000 153.0333
## coords.x2  -37.81667  57.0500
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 54

Los objetos SpatialPoints pueden tener más de dos dimensiones, pero los métodos de visualización usan solo las dos primeras.

3.4.1 Métodos

Hay métodos disponibles para acceder a los valores de las ranuras de los objetos espaciales.

  • bbox: devuelve los límites del objeto y se utiliza tanto para preparar métodos de visualización como internamente en el manejo de objetos de datos.

Si queremos tomar un subconjunto de los puntos en un objeto SpatialPoints, entonces bbox se restablece.

sp::bbox(CRAN_sp)
##                  min      max
## coords.x1 -122.95000 153.0333
## coords.x2  -37.81667  57.0500
  • proj4string: reporta la proyección contenida como un objeto CRS en el slot proj4string del objeto, también le permite al usuario alterar el valor actual:
sp::proj4string(CRAN_sp)
## [1] "+proj=longlat +ellps=WGS84 +no_defs"
sp::proj4string(CRAN_sp) <- llCRS
  • coordinates: extrae las coordenadas del objeto SpatialPoints.
brazil <- which(CRAN_df$loc == "Brazil")
brazil
## [1] 4 5 6 7 8
sp::coordinates(CRAN_sp)[brazil, ]
##   coords.x1 coords.x2
## 4 -49.26667 -25.41667
## 5 -42.86667 -20.75000
## 6 -43.20000 -22.90000
## 7 -47.63333 -22.71667
## 8 -46.63333 -23.53333

También se puede hace indexación usando corechetes:

sp::summary(CRAN_sp[brazil, ])
## Object of class SpatialPoints
## Coordinates:
##                 min       max
## coords.x1 -49.26667 -42.86667
## coords.x2 -25.41667 -20.75000
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 5
south_of_equator <- which(sp::coordinates(CRAN_sp)[, 2] < 0)
sp::summary(CRAN_sp[-south_of_equator, ])
## Object of class SpatialPoints
## Coordinates:
##               min    max
## coords.x1 -122.95 140.10
## coords.x2   24.15  57.05
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 45
  • summary: reporta el número de entidades espaciales, información de la proyección y los límites.

  • print: da un vistazo del objeto.

3.4.2 Data Frames para puntos espaciales

Los SpatialPoints se pueden comportar como un data.frame. Se usan números en secuencia para indexar los puntos y filas en nuestro data frame.

str(row.names(CRAN_df))
##  chr [1:54] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" ...

Si la matriz de coordenadas de puntos tiene nombres de fila y el argumento match.ID se establece en su valor predeterminado de TRUE, entonces los nombres de fila de la matriz se verifican con los nombres de fila del data frame:

  • Si coinciden, pero no están en el mismo orden, las filas del data.frame se reordenan para adaptarse a los puntos.
  • Si no coinciden, no se construye ningún SpatialPointsDataFrame.
CRAN_spdf1 <- sp::SpatialPointsDataFrame(CRAN_mat, CRAN_df, proj4string = llCRS, match.ID = TRUE)
CRAN_spdf1[4, ]
## class       : SpatialPointsDataFrame 
## features    : 1 
## extent      : -49.26667, -49.26667, -25.41667, -25.41667  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +ellps=WGS84 +no_defs 
## variables   : 6
## names       :    place,   north,    east,    loc,              long,               lat 
## value       : Curitiba, 25d25'S, 49d16'W, Brazil, -49.2666666666667, -25.4166666666667

Podemos acceder a los datos como si usáramos data.frames.

str(CRAN_spdf1$loc)
##  chr [1:54] "Australia" "Australia" "Austria" "Brazil" "Brazil" "Brazil" "Brazil" "Brazil" ...
str(CRAN_spdf1[["loc"]])
##  chr [1:54] "Australia" "Australia" "Austria" "Brazil" "Brazil" "Brazil" "Brazil" "Brazil" ...

Pero si tenemos valores de ID que no coinciden, creados al unir pares de letras y muestrear un número apropiado de ellas, el resultado es un error:

CRAN_df1 <- CRAN_df
row.names(CRAN_df1) <- sample(c(outer(letters, letters,paste, sep = "")), nrow(CRAN_df1))
CRAN_spdf3 <- sp::SpatialPointsDataFrame(CRAN_mat, CRAN_df1, proj4string = llCRS, match.ID = TRUE)

Un objeto SpatialPointsDataFrame hereda información contenida en el objeto Spatial:

getClass("SpatialPointsDataFrame")
## Class "SpatialPointsDataFrame" [package "sp"]
## 
## Slots:
##                                                                   
## Name:         data  coords.nrs      coords        bbox proj4string
## Class:  data.frame     numeric      matrix      matrix         CRS
## 
## Extends: 
## Class "SpatialPoints", directly
## Class "Spatial", by class "SpatialPoints", distance 2
## Class "SpatialVector", by class "SpatialPoints", distance 2
## Class "SpatialPointsNULL", by class "SpatialPoints", distance 2
## 
## Known Subclasses: 
## Class "SpatialPixelsDataFrame", directly, with explicit coerce

Otra forma de crear un SpatialPointsDataFrame:

CRAN_spdf4 <- sp::SpatialPointsDataFrame(CRAN_sp, CRAN_df)
all.equal(CRAN_spdf4, CRAN_spdf1)
## [1] TRUE

La función de asignación de coordenadas puede tomar una matriz de coordenadas con el mismo número de filas que el data frame en el lado derecho:

CRAN_df0 <- CRAN_df
sp::coordinates(CRAN_df0) <- CRAN_mat
sp::proj4string(CRAN_df0) <- llCRS
all.equal(CRAN_df0, CRAN_spdf1)
## [1] TRUE
str(CRAN_df0, max.level = 2)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  54 obs. of  6 variables:
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:54, 1:2] 153 145 16.3 -49.3 -42.9 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ bbox       : num [1:2, 1:2] -123 -37.8 153 57
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

Números en secuencia se pueden ingresar en el data frame para que sea posible rastrear los puntos en orden, por ejemplo, como parte de un objeto SpatialLines. Por ejemplo:

uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/seamap105_mod.csv"

turtle_df <- read.csv(url(uu))
summary(turtle_df)
##        id              lat             lon            obs_date        
##  Min.   :  1.00   Min.   :21.57   Min.   :-179.88   Length:394        
##  1st Qu.: 99.25   1st Qu.:24.36   1st Qu.:-147.38   Class :character  
##  Median :197.50   Median :25.64   Median :-119.64   Mode  :character  
##  Mean   :197.50   Mean   :27.21   Mean   : -21.52                     
##  3rd Qu.:295.75   3rd Qu.:27.41   3rd Qu.: 153.66                     
##  Max.   :394.00   Max.   :39.84   Max.   : 179.93

Antes de crear un SpatialPointsDataFrame, marcaremos el tiempo de las observaciones y reordenaremos el data frame de entrada por el tiempo para que sea más fácil agregar meses.

timestamp <- as.POSIXlt(strptime(as.character(turtle_df$obs_date), "%m/%d/%Y %H:%M:%S"), "GMT")
turtle_df1 <- data.frame(turtle_df, timestamp = timestamp)
turtle_df1$lon <- ifelse(turtle_df1$lon < 0, turtle_df1$lon + 360, turtle_df1$lon)
turtle_sp <- turtle_df1[order(turtle_df1$timestamp),]
sp::coordinates(turtle_sp) <- c("lon", "lat")
sp::proj4string(turtle_sp) <- sp::CRS("+proj=longlat +ellps=WGS84")

Ahora vemos el mapa:

# library(maptools)
gshhs.c.b <- system.file("share/gshhs_c.b", package="maptools")
pac <- maptools::Rgshhs(gshhs.c.b, level=1, xlim=c(130,250), ylim=c(15,60), verbose=FALSE)

sp::plot(pac$SP, axes=TRUE, col="khaki2", xaxs="i", yaxs="i")
sp::plot(turtle_sp, add=TRUE)
m_rle <- rle(months(turtle_sp$timestamp))
clen <- cumsum(m_rle$lengths[-length(m_rle$lengths)])-1
crds <- sp::coordinates(turtle_sp)
text(crds[clen,], labels=m_rle$values[-1], pos=3, offset=1.5, srt=45)

3.5 Líneas espaciales

Un objeto Línea es una matriz de coordenadas 2D, sin valores NA

un objeto geométrico unidimensional normalmente almacenado como una secuencia de puntos. . . con interpolación lineal entre puntos (Herring and others (2011))

La interpolación lineal significa que asumimos que las coordenadas intermedias no observadas pueden interpolarse usando una línea recta entre las coordenadas en cualquier extremo del segmento de línea.

Una lista de objetos de Line forma el slot de líneas de un objeto de Lines

También se requiere una etiqueta de identificación, y se utilizará para construir objetos SpatialLines.

getClass("Line")
## Class "Line" [package "sp"]
## 
## Slots:
##              
## Name:  coords
## Class: matrix
## 
## Known Subclasses: "Polygon"
getClass("Lines")
## Class "Lines" [package "sp"]
## 
## Slots:
##                           
## Name:      Lines        ID
## Class:      list character

Objeto SpatialLines que contiene el cuadro delimitador y la información de proyección para la lista de objetos Lines almacenados en su slot de líneas.

getClass("SpatialLines")
## Class "SpatialLines" [package "sp"]
## 
## Slots:
##                                           
## Name:        lines        bbox proj4string
## Class:        list      matrix         CRS
## 
## Extends: "Spatial", "SpatialVector", "SpatialLinesNULL"
## 
## Known Subclasses: "SpatialLinesDataFrame"

Examinemos un ejemplo creado a partir de líneas de la base de datos mundial del paquete de maps, y convertido a un objeto SpatialLines usando la función map2SpatialLines en maptools.

# library(maps)
japan <- maps::map("world", "japan", plot = FALSE)
p4s <- sp::CRS("+proj=longlat +ellps=WGS84")
# library(maptools)
SLjapan <- maptools::map2SpatialLines(japan, proj4string = p4s)
str(SLjapan, max.level = 2)
## Formal class 'SpatialLines' [package "sp"] with 3 slots
##   ..@ lines      :List of 34
##   ..@ bbox       : num [1:2, 1:2] 123.7 24.3 145.8 45.5
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

Los objetos SpatialLines y SpatialPolygons son muy similares, se suele usar lapply o sapply para explorar los contenidos de estos objetos en combinación con slot

Lines_len <- sapply(slot(SLjapan, "lines"), function(x) length(slot(x, "Lines")))
table(Lines_len)
## Lines_len
##  1 
## 34

Podemos usar la función ContourLines2SLDF incluida en maptools en nuestro siguiente ejemplo, convirtiendo los datos devueltos por la función gráfica base contourLines en un objeto SpatialLinesDataFrame

volcano_sl <- maptools::ContourLines2SLDF(contourLines(z = volcano))
t(slot(volcano_sl, "data"))
##       C_1   C_2   C_3   C_4   C_5   C_6   C_7   C_8   C_9   C_10 
## level "100" "110" "120" "130" "140" "150" "160" "170" "180" "190"
sp::plot(volcano_sl)

Podemos ver que hay diez etiquetas de nivel de contorno separadas en el slot, almacenadas como un factor en el data frame del slot de datos del objeto.

llCRS <- sp::CRS("+proj=longlat +ellps=WGS84")
uu <- ("https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/auckland_mapgen.dat")
auck_shore <- maptools::MapGen2SL(uu, llCRS)
sp::summary(auck_shore)
## Object of class SpatialLines
## Coordinates:
##     min   max
## x 174.2 175.3
## y -37.5 -36.5
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]

3.6 Polígonos espaciales

Es una secuencia de coordenadas de puntos donde el primer punto es el mismo que el último punto.

Para tener un conjunto de datos para usar en polígonos, primero identificamos las líneas importadas arriba que representan la costa alrededor de Auckland. Muchas son islas, por lo que tienen una primera y una última coordenadas idénticas.

lns <- slot(auck_shore, "lines")
table(sapply(lns, function(x) length(slot(x, "Lines"))))
## 
##  1 
## 80
islands_auck <- sapply(lns, function(x) {
  crds <- slot(slot(x, "Lines")[[1]], "coords")
  identical(crds[1,], crds[nrow(crds),]) # miramos que coincida la primera y última coordenada
})
table(islands_auck)
## islands_auck
## FALSE  TRUE 
##    16    64

Veamos el objeto polígono:

getClass("Polygon")
## Class "Polygon" [package "sp"]
## 
## Slots:
##                                               
## Name:    labpt    area    hole ringDir  coords
## Class: numeric numeric logical integer  matrix
## 
## Extends: "Line"

La clase Polygon amplía la clase Line agregando los slots necesarios para los polígonos y verificando que la primera y la última coordenadas sean idénticas.

  • labpt: centroide del polígono
  • area: área del polígono
  • hole: si el polígono es declarado con un agujero o no.
  • ringDir: dirección del anillo.

hole y ringDir se incluyen en los objetos Polygon como heurísticas para abordar algunas de las dificultades que surgen de que R no sea un GIS. El enfoque que se ha elegido en sp es utilizar dos marcadores que se encuentran comúnmente en la práctica, marcar polígonos como agujeros con una bandera lógica (VERDADERO / FALSO), el hole, y usar la dirección del anillo: los anillos en el sentido de las agujas del reloj se consideran que no son agujeros, en sentido antihorario como agujeros. Esto es necesario porque la representación no topológica de polígonos no tiene una manera fácil de saber que un polígono representa un límite interno de un polígono circundante, un agujero o un lago.

Veamos el conjunto de polígonos:

getClass("Polygons")
## Class "Polygons" [package "sp"]
## 
## Slots:
##                                                         
## Name:   Polygons plotOrder     labpt        ID      area
## Class:      list   integer   numeric character   numeric

Veamos ahora los polígonos espaciales:

getClass("SpatialPolygons")
## Class "SpatialPolygons" [package "sp"]
## 
## Slots:
##                                                       
## Name:     polygons   plotOrder        bbox proj4string
## Class:        list     integer      matrix         CRS
## 
## Extends: "Spatial", "SpatialVector", "SpatialPolygonsNULL"
## 
## Known Subclasses: 
## Class "SpatialPolygonsDataFrame", directly, with explicit coerce

Es un conjunto de polígonos donde ya podemos tener la información de bbox y de CRS. Los polígonos se almacenan del área más grande a la más pequeña.

Al elegir solo las líneas en el conjunto de datos de la costa de Auckland que son polígonos cerrados, podemos construir un objeto SpatialPolygons.

islands_sl <- auck_shore[islands_auck]
list_of_Lines <- slot(islands_sl, "lines")
islands_sp <- sp::SpatialPolygons(lapply(list_of_Lines, function(x) {
    sp::Polygons(list(sp::Polygon(slot(slot(x, "Lines")[[1]], "coords"))),
      ID=slot(x, "ID"))
  }),
  proj4string=sp::CRS("+proj=longlat +ellps=WGS84"))
sp::summary(islands_sp)
## Object of class SpatialPolygons
## Coordinates:
##         min       max
## x 174.30297 175.22791
## y -37.43877 -36.50033
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]

Veamos el orden en el que se grafica:

slot(islands_sp, "plotOrder")
##  [1] 45 54 37 28 38 27 12 11 59 53  5 25 26 46  7 55 17 34 30 16  6 43 14 40 32 19 61 42 15 50 21 18 62
## [34] 23 22 29 24 44 13  2 36  9 63 58 56 64 52 39 51  1  8  3  4 20 47 35 41 48 60 31 49 57 10 33

Verificamos que se grafican por tamaño del área:

order(sapply(slot(islands_sp, "polygons"), function(x) slot(x, "area")), decreasing = TRUE)
##  [1] 45 54 37 28 38 27 12 11 59 53  5 25 26 46  7 55 17 34 30 16  6 43 14 40 32 19 61 42 15 50 21 18 62
## [34] 23 22 29 24 44 13  2 36  9 63 58 56 64 52 39 51  1  8  3  4 20 47 35 41 48 60 31 49 57 10 33
sp::plot(auck_shore)
legend("bottomleft", legend="a)", bty="n")

sp::plot(auck_shore)
sp::plot(islands_sp, add=TRUE, col="grey")
legend("bottomleft", legend="b)", bty="n")

3.6.1 Data Frames para polígonos

Los objetos SpatialPolygonsDataFrame reúnen las representaciones espaciales de los polígonos con datos.

Las etiquetas de identificación de los polígonos en el slot de polígono de un objeto SpatialPolygons coinciden con los nombres de fila del data.frame para asegurarse de que las filas de datos correctas estén asociadas con los objetos espaciales correctos.

Como ejemplo, tomamos un conjunto de puntajes por el estado de EE. UU. De 1999 Scholastic Aptitude Test (SAT)

# library(maps)
state.map <- maps::map("state", plot = FALSE, fill = TRUE)
IDs <- sapply(strsplit(state.map$names, ":"), function(x) x[1])
# library(maptools)
state.sp <- maptools::map2SpatialPolygons(state.map, IDs = IDs,proj4string = sp::CRS("+proj=longlat +ellps=WGS84"))

Luego, podemos usar la coincidencia de etiquetas de identificación para unir por las filas del data.frame a los SpatialPolygons. Aquí, hacemos un subconjunto de las filas coincidentes del data.frame, para asegurarnos de que una fila corresponda a cada objeto Polygons, para lograr una coincidencia uno a uno:

uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/state.sat.data_mod.txt"
sat <- read.table(uu, row.names = 1, header = TRUE)
str(sat)
## 'data.frame':    52 obs. of  4 variables:
##  $ oname : chr  "ala" "alaska" "ariz" "ark" ...
##  $ vscore: int  561 516 524 563 497 536 510 503 494 499 ...
##  $ mscore: int  555 514 525 556 514 540 509 497 478 498 ...
##  $ pc    : int  9 50 34 6 49 32 80 67 77 53 ...

Creamos la variable de match:

id <- match(row.names(sat), row.names(state.sp))
row.names(sat)[is.na(id)]
## [1] "alaska" "hawaii" "usa"

Agregamos la variable:

sat1 <- sat[!is.na(id), ]
state.spdf <- sp::SpatialPolygonsDataFrame(state.sp, sat1)
str(slot(state.spdf, "data"))
## 'data.frame':    49 obs. of  4 variables:
##  $ oname : chr  "ala" "ariz" "ark" "calif" ...
##  $ vscore: int  561 524 563 497 536 510 503 494 499 487 ...
##  $ mscore: int  555 525 556 514 540 509 497 478 498 482 ...
##  $ pc    : int  9 34 6 49 32 80 67 77 53 63 ...
str(state.spdf, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  49 obs. of  4 variables:
##   ..@ polygons   :List of 49
##   ..@ plotOrder  : int [1:49] 42 25 4 30 27 2 5 49 36 22 ...
##   ..@ bbox       : num [1:2, 1:2] -124.7 25.1 -67 49.4
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

Para quitar un distrito:

DC <- "district of columbia"
not_dc <- !(row.names(state.spdf) == DC)
state.spdf1 <- state.spdf[not_dc, ]
dim(state.spdf1)
## [1] 48  4
summary(state.spdf1)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##          min       max
## x -124.68134 -67.00742
## y   25.12993  49.38323
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Data attributes:
##     oname               vscore          mscore            pc       
##  Length:48          Min.   :479.0   Min.   :475.0   Min.   : 4.00  
##  Class :character   1st Qu.:506.2   1st Qu.:505.2   1st Qu.: 9.00  
##  Mode  :character   Median :530.5   Median :532.0   Median :28.50  
##                     Mean   :534.6   Mean   :534.9   Mean   :35.58  
##                     3rd Qu.:563.0   3rd Qu.:558.5   3rd Qu.:63.50  
##                     Max.   :594.0   Max.   :605.0   Max.   :80.00

3.7 Objetos grid

  • Las cuadrículas son objetos regulares que requieren mucha menos información para definir su estructura.

  • Una vez que se conoce el único punto de origen, la extensión de la cuadrícula puede estar dada por la resolución de la celda y el número de filas y columnas presentes en la cuadrícula completa.

getClass("GridTopology")
## Class "GridTopology" [package "sp"]
## 
## Slots:
##                                                             
## Name:  cellcentre.offset          cellsize         cells.dim
## Class:           numeric           numeric           integer

Como ejemplo, creamos un objeto GridTopology a partir del bbox del conjunto de datos vectoriales de la isla Manitoulin. Si elegimos un tamaño de celda de 0.01º en cada dirección, podemos compensar el centro de celda suroeste para asegurarnos de que al menos toda el área esté cubierta y encontrar un número adecuado de celdas en cada dimensión.

uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/high.RData"
load(url(uu))
manitoulin_sp <- high$SP
bb <- sp::bbox(manitoulin_sp)
bb
##      min   max
## r1 277.0 278.0
## r2  45.7  46.2
cs <- c(0.01, 0.01)
cc <- bb[, 1] + (cs/2)
cd <- ceiling(diff(t(bb))/cs)
manitoulin_grd <- sp::GridTopology(cellcentre.offset = cc,cellsize = cs, cells.dim = cd)
manitoulin_grd
##                        r1     r2
## cellcentre.offset 277.005 45.705
## cellsize            0.010  0.010
## cells.dim         100.000 50.000

El objeto describe la cuadrícula por completo y se puede usar para construir un objeto SpatialGrid. Un objeto SpatialGrid contiene objetos GridTopology y Spatial.

getClass("SpatialGrid")
## Class "SpatialGrid" [package "sp"]
## 
## Slots:
##                                              
## Name:          grid         bbox  proj4string
## Class: GridTopology       matrix          CRS
## 
## Extends: "Spatial"
## 
## Known Subclasses: "SpatialGridDataFrame"

Usando el objeto GridTopology creado anteriormente, y pasando por el sistema de referencia de coordenadas de los datos GSHHS originales, el bbox se crea automáticamente, como vemos en el resumen del objeto.

p4s <- sp::CRS(sp::proj4string(manitoulin_sp))
manitoulin_SG <- sp::SpatialGrid(manitoulin_grd, proj4string = p4s)
sp::summary(manitoulin_SG)
## Object of class SpatialGrid
## Coordinates:
##      min   max
## r1 277.0 278.0
## r2  45.7  46.2
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Grid attributes:
##    cellcentre.offset cellsize cells.dim
## r1           277.005     0.01       100
## r2            45.705     0.01        50

3.8 Objetos raster

Los ráster con datos están representadas por las coordenadas del centro de la celda y por el número de secuencia de la celda.

Un avance del paquete raster es que los objetos ráster pueden mantenerse en el disco en lugar de en la memoria.

# library(raster)
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/70042108.tif"
r <- raster::raster(uu)
class(r)
## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"
raster::inMemory(r)
## [1] FALSE
object.size(r)
## 13536 bytes
raster::cellStats(r, max)
## [1] 686
raster::cellStats(r, min)
## [1] -3.402823e+38

Necesitamos eliminar los valores menores o iguales a cero, ya que las elevaciones solo se registran en tierra.

out <- raster::raster(r)
bs <- raster::blockSize(out)
out <- raster::writeStart(out, filename = tempfile(), overwrite = TRUE)
for (i in 1:bs$n) {
v <- raster::getValues(r, row = bs$row[i], nrows = bs$nrows[i])
v[v <= 0] <- NA
raster::writeValues(out, v, bs$row[i])
}
out <- raster::writeStop(out)
raster::cellStats(out, min)
## [1] 1
raster::cellStats(out, max)
## [1] 686
sp::plot(out, col = terrain.colors(100))

El RasterLayer puede convertirse en SpatialGridDataFrame y al revés:

r1 <- as(out, "SpatialGridDataFrame")
sp::summary(r1)
## Object of class SpatialGridDataFrame
## Coordinates:
##      min   max
## s1 174.2 175.3
## s2 -37.5 -36.5
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Grid attributes:
##    cellcentre.offset     cellsize cells.dim
## s1         174.20042 0.0008333333      1320
## s2         -37.49958 0.0008333333      1200
## Data attributes:
##      layer       
##  Min.   :  1.0   
##  1st Qu.: 23.0   
##  Median : 53.0   
##  Mean   : 78.1   
##  3rd Qu.:106.0   
##  Max.   :686.0   
##  NA's   :791938
r2 <- as(r1, "RasterLayer")
sp::summary(r2)
##          layer
## Min.         1
## 1st Qu.     23
## Median      53
## 3rd Qu.    106
## Max.       686
## NA's    791938

La formalidad que exige el tratar adecuadamente los objetos espaciales toma relevancia al momento de generar visualizaciones.

3.9 Referencias