3 Objetos espaciales

Tomado de Pebesma and Bivand (2023), espefíficamente de las viñetas(https://r-spatial.github.io/sf/index.html)

Librerías usadas:

• sf

3.1 Atributos simples

Atributos simples o acceso simple a atributos se refiere a un estándar formal (ISO 19125-1:2004) que describe:

1. cómo los objetos del mundo real se pueden representar en computadoras, con énfasis en la geometría espacial de estos objetos.
2. cómo dichos objetos pueden almacenarse y extraerse de bases de datos, y
3. qué operaciones geométricas deben definirse para ellos.

El estándar se implementa ampliamente en bases de datos espaciales (como PostGIS), Sistemas de Información Geográfica -SIG comerciales (por ejemplo, ESRI ArcGIS) y forma la base de datos vectoriales para bibliotecas como GDAL. Un subconjunto de funciones simples forma el estándar GeoJSON.

R tiene clases bien soportadas para almacenar datos espaciales (sp) y conectarse a los entornos mencionados anteriormente (rgdal,rgeos), pero hasta ahora ha carecido de una implementación completa de características simples, lo que hace que las conversiones a veces sean complicadas, ineficientes o incompletas. El paquete sf intenta llenar este vacío y apunta a suceder al sp a largo plazo.

3.2 ¿Qué es una característica (feature)?

Una característica se considera una cosa o un objeto del mundo real, como un edificio o un árbol. Como ocurre con los objetos, a menudo se componen de otros objetos. Este también es el caso de las características: un conjunto de características puede formar una sola característica. Una masa forestal puede ser una característica, un bosque puede ser una característica, una ciudad puede ser una característica. Un píxel de una imagen de satélite puede ser una característica, y una imagen completa también puede ser una característica.

Las entidades tienen una geometría que describe en qué lugar de la Tierra se encuentra la entidad y tienen atributos que describen otras propiedades. La geometría de un árbol puede ser la delimitación de su copa, de su tallo o el punto que indica su centro. Otras propiedades pueden incluir su altura, color, diámetro a la altura del pecho en una fecha determinada, etc.

3.2.1 Dimensiones

Todas las geometrías están compuestas de puntos. Los puntos son coordenadas en un espacio de 2, 3 o 4 dimensiones. Todos los puntos de una geometría tienen la misma dimensionalidad. Además de las coordenadas \(X\) e \(Y\), existen dos dimensiones adicionales opcionales:

• una coordenada \(Z\) que denota altura.
• una coordenada \(M\) (rara vez utilizada), que denota alguna medida asociada con el punto, en lugar de con la característica en su conjunto; algunos ejemplos podrían ser el tiempo de medición o el error de medición de las coordenadas.

Los cuatro casos posibles se denota como:

1. Los puntos bidimensionales se refieren a x e y, este y norte, o longitud y latitud, nos referimos a ellos como \(XY\).
2. puntos tridimensionales como \(XYZ\).
3. puntos tridimensionales como \(XYM\).
4. puntos en cuatro dimensiones como \(XYZM\).

3.2.2 Tipos de geometría de entidades simples

Los siguientes siete tipos de funciones simples son los más comunes y, por ejemplo, los únicos que se utilizan para GeoJSON:

type	description
POINT	geometría de dimensión cero que contiene un solo punto
LINESTRING	secuencia de puntos conectados por líneas rectas que no se cruzan entre sí; geometría unidimensional
POLYGON	geometría con área positiva (bidimensional); la secuencia de puntos forma un anillo cerrado que no se cruza entre sí; el primer anillo indica el anillo exterior, cero o más anillos posteriores indican agujeros en este anillo exterior
MULTIPOINT	set of points; a MULTIPOINT is simple if no two Points in the MULTIPOINT are equal
MULTILINESTRING	set of linestrings
MULTIPOLYGON	set of polygons
GEOMETRYCOLLECTION	set of geometries of any type except GEOMETRYCOLLECTION

Cada uno de los tipos de geometría también puede ser un conjunto vacío, que contiene coordenadas cero (para PUNTO, el estándar no está claro cómo representar la geometría vacía). Se puede considerar que las geometrías vacías son análogas a los atributos faltantes (NA), los valores NULL o las listas vacías.

3.2.3 Sistema de Referencia de Coordenadas

Coordenadas polares y cartesianas:

\[x = r \cos \phi, y = r \sin \phi\] \[r = \sqrt{x^2+y^2}, \phi = atan2(x,y)\]

Donde atan2: La función arcotangente de dos parámetros devuelve el ángulo formado entre el eje \(x\) positivo y la recta que conecta el origen con un punto de coordenadas \((x,y)\neq(0,0)\) del plano euclidiano, expresado en radianes.

Coordenadas esféricas:

En tres dimensiones las coordenadas cartesianas se expresan como \((x,y,z)\), las coordenadas esféricas se expresan con \((r,\lambda,\phi)\) donde

• \(r\) es el radio de la esfera,
• \(\lambda\) es la longitud, medida en el plano \((x,y)\) en contra de las manecillas del reloj desde el eje \(x\) positivo
• \(\phi\) es la latitud, el ángulo entre el vector y el plano \((x,y)\).

• \(-180º\leq\lambda\leq180º\) (o \(0º\leq\lambda\leq360º\)).
• \(-90º\leq\phi\leq90º\).
• Si estamos interesados en un punto de la esfera, se suele omitir \(r\) porque en ese caso \((\lambda,\phi)\) son suficientes para identificar el punto.

Coordenadas en la elipse:

Hay dos maneras en las que un ángulo puede expresarse:

• medido desde el centro de la elipse (\(\psi\)) o,
• medido perpendicular a la tangente de la elipse en el punto de medición (\(\phi\)).

• Un sistema de coordenadas es un conjunto de reglas matemáticas para especificar cómo se asignarán las coordenadas a los puntos.
• Un datum es un parámetro o conjunto de parámetros que definen la posición del origen, la escala, y la orientación del sistema de coordenadas.
• Un geodetic datum o datum geodésico es un datum que describe la relación de un sistema de coordenadas a la tierra.
• Un sistema de referencia de coordenadas es un sistema de coordenadas que está relacionado a un objeto (a la tierra en el caso geodésico) por un datum.

Video: https://www.youtube.com/watch?v=SaFrbKVC1Yw

Las coordenadas sólo se pueden colocar en la superficie de la Tierra cuando se conoce su sistema de referencia de coordenadas (CRS); este puede ser

• un CRS esferoide como WGS84,
• un CRS bidimensional (cartesiano) proyectado como una zona UTM o Web Mercator,
• o un CRS en tres dimensiones, o incluyendo el tiempo.

3.3 Cómo se organizan los atributos simples

El paquete sf representa características simples como objetos R nativos. De manera similar a PostGIS, todas las funciones y métodos en sf que operan con datos espaciales tienen el prefijo st_, que se refiere al tipo espacial; esto hace que sean fáciles de encontrar.

Las funciones simples se implementan como datos nativos de R, utilizando estructuras de datos simples (clases S3, listas, matrices, vectores). El uso típico implica leer, manipular y escribir conjuntos de características, con atributos y geometrías.

Como los atributos generalmente se almacenan en objetos data.frame (o el muy similar tbl_df), también almacenaremos geometrías de características en una columna data.frame.

Dado que las geometrías no tienen un solo valor, se colocan en una columna de lista, una lista de longitud igual al número de registros en el marco de datos, y cada elemento de la lista contiene la geometría de característica simple de esa característica.

Las tres clases utilizadas para representar características simples son:

• sf, la tabla (data.frame) con atributos de características y geometrías de características, que contiene
• sfc, la lista-columna con las geometrías para cada característica (registro), que se compone de
• sfg, la geometría de característica de una característica simple individual.

library(sf)
nc <- st_read(system.file("shape/nc.shp", package="sf"))

## Reading layer `nc' from data source 
##   `/Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/sf/shape/nc.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27

Veamos su clase

class(nc)

## [1] "sf"         "data.frame"

Para identificar el nombre de la columna que tiene las geometrías:

attr(nc, "sf_column")

## [1] "geometry"

Si se imprime tres valores

print(nc[9:15], n = 3)

En la salida vemos:

• en verde una característica simple: un único registro, o fila del data.frame, que consta de atributos y geometría
• en azul una característica simple (un objeto de clase sfg)
• en rojo, una columna de lista de características simples (un objeto de clase sfc, que es una columna en el data.frame)

Los métodos para objetos sf son:

methods(class = "sf")

##   [1] [                            [[<-                        
##   [3] [<-                          $<-                         
##   [5] aggregate                    anti_join                   
##   [7] arrange                      as.data.frame               
##   [9] cbind                        coerce                      
##  [11] crs                          dbDataType                  
##  [13] dbWriteTable                 distance                    
##  [15] distinct                     dplyr_reconstruct           
##  [17] duplicated                   extent                      
##  [19] extract                      filter                      
##  [21] full_join                    group_by                    
##  [23] group_split                  identify                    
##  [25] idw                          initialize                  
##  [27] inner_join                   krige                       
##  [29] krige.cv                     left_join                   
##  [31] lines                        mask                        
##  [33] merge                        mutate                      
##  [35] plot                         print                       
##  [37] raster                       rasterize                   
##  [39] rbind                        rename_with                 
##  [41] rename                       right_join                  
##  [43] rowwise                      sample_frac                 
##  [45] sample_n                     select                      
##  [47] semi_join                    show                        
##  [49] slice                        slotsFromS3                 
##  [51] st_agr                       st_agr<-                    
##  [53] st_area                      st_as_s2                    
##  [55] st_as_sf                     st_as_sfc                   
##  [57] st_bbox                      st_boundary                 
##  [59] st_break_antimeridian        st_buffer                   
##  [61] st_cast                      st_centroid                 
##  [63] st_collection_extract        st_concave_hull             
##  [65] st_convex_hull               st_coordinates              
##  [67] st_crop                      st_crs                      
##  [69] st_crs<-                     st_difference               
##  [71] st_drop_geometry             st_filter                   
##  [73] st_geometry                  st_geometry<-               
##  [75] st_inscribed_circle          st_interpolate_aw           
##  [77] st_intersection              st_intersects               
##  [79] st_is_valid                  st_is                       
##  [81] st_join                      st_line_merge               
##  [83] st_m_range                   st_make_valid               
##  [85] st_minimum_rotated_rectangle st_nearest_points           
##  [87] st_node                      st_normalize                
##  [89] st_point_on_surface          st_polygonize               
##  [91] st_precision                 st_reverse                  
##  [93] st_sample                    st_segmentize               
##  [95] st_set_precision             st_shift_longitude          
##  [97] st_simplify                  st_snap                     
##  [99] st_sym_difference            st_transform                
## [101] st_triangulate_constrained   st_triangulate              
## [103] st_union                     st_voronoi                  
## [105] st_wrap_dateline             st_write                    
## [107] st_z_range                   st_zm                       
## [109] summarise                    transform                   
## [111] transmute                    ungroup                     
## see '?methods' for accessing help and source code

Podemos transformar un objeto sf a data.frame (se pierde las características sf por lo que ya no aplicarían los métodos anteriores)

nc.no_sf <- as.data.frame(nc)
class(nc.no_sf)

## [1] "data.frame"

Cheatsheet: https://github.com/rstudio/cheatsheets/blob/main/sf.pdf

3.4 sfc: la columna con la geometría simple

La columna en el sf data.frame que contiene las geometrías es una lista, de clase sfc. Podemos recuperar la columna de lista de geometría en este caso mediante nc$geom o nc[[15]], pero la forma más general usa `st_geometry()```:

nc$geometry

## Geometry set for 100 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## First 5 geometries:

## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...

## MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...

## MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...

## MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...

## MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...

(nc_geom <- st_geometry(nc))

## Geometry set for 100 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## First 5 geometries:

## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...

## MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...

## MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...

## MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...

## MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...

Las geometrías se imprimen en forma abreviada, pero podemos ver una geometría completa seleccionándola, p.e. el primero por:

nc_geom[[1]]

## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.34069, -81.74107 36.39178, -81.69828 36.47178, -81.7028 36.51934, -81.67 36.58965, -81.3453 36.57286, -81.34754 36.53791, -81.32478 36.51368, -81.31332 36.4807, -81.26624 36.43721, -81.26284 36.40504, -81.24069 36.37942, -81.23989 36.36536, -81.26424 36.35241, -81.32899 36.3635, -81.36137 36.35316, -81.36569 36.33905, -81.35413 36.29972, -81.36745 36.2787, -81.40639 36.28505, -81.41233 36.26729, -81.43104 36.26072, -81.45289 36.23959, -81.47276 36.23436)))

La clumna de geometrías tiene su propia clase:

class(nc_geom)

## [1] "sfc_MULTIPOLYGON" "sfc"

Con los siguientes métodos:

methods(class = 'sfc')

##  [1] [                            [<-                         
##  [3] as.data.frame                c                           
##  [5] coerce                       format                      
##  [7] identify                     initialize                  
##  [9] Ops                          print                       
## [11] rep                          show                        
## [13] slotsFromS3                  st_area                     
## [15] st_as_binary                 st_as_grob                  
## [17] st_as_s2                     st_as_sf                    
## [19] st_as_text                   st_bbox                     
## [21] st_boundary                  st_break_antimeridian       
## [23] st_buffer                    st_cast                     
## [25] st_centroid                  st_collection_extract       
## [27] st_concave_hull              st_convex_hull              
## [29] st_coordinates               st_crop                     
## [31] st_crs                       st_crs<-                    
## [33] st_difference                st_geometry                 
## [35] st_inscribed_circle          st_intersection             
## [37] st_intersects                st_is_valid                 
## [39] st_is                        st_line_merge               
## [41] st_m_range                   st_make_valid               
## [43] st_minimum_rotated_rectangle st_nearest_points           
## [45] st_node                      st_normalize                
## [47] st_point_on_surface          st_polygonize               
## [49] st_precision                 st_reverse                  
## [51] st_sample                    st_segmentize               
## [53] st_set_precision             st_shift_longitude          
## [55] st_simplify                  st_snap                     
## [57] st_sym_difference            st_transform                
## [59] st_triangulate_constrained   st_triangulate              
## [61] st_union                     st_voronoi                  
## [63] st_wrap_dateline             st_write                    
## [65] st_z_range                   st_zm                       
## [67] str                          summary                     
## [69] type_sum                     vec_cast.sfc                
## [71] vec_ptype2.sfc              
## see '?methods' for accessing help and source code

Atributos

attributes(nc_geom)

## $n_empty
## [1] 0
## 
## $crs
## Coordinate Reference System:
##   User input: NAD27 
##   wkt:
## GEOGCRS["NAD27",
##     DATUM["North American Datum 1927",
##         ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4267]]
## 
## $class
## [1] "sfc_MULTIPOLYGON" "sfc"             
## 
## $precision
## [1] 0
## 
## $bbox
##      xmin      ymin      xmax      ymax 
## -84.32385  33.88199 -75.45698  36.58965

3.5 sfg Geometría de caracterísitia simple

Los objetos de geometría de característica simple (sfg) llevan la geometría de una sola característica: un punto, cadena lineal (LINESTRING) o polígono.

Las características simples se implementan en R con tipos nativos, siguiendo las siguientes reglas:

1. un solo POINT es un vector numérico
2. un conjunto de puntos (como LINESTRING o POLYGON) es una matriz, donde cada fila contiene un punto.
3. cualquier otro conjunto es list.

Como ejemplo:

(x <- st_point(c(1,2)))

## POINT (1 2)

(x <- st_point(c(1,2,3)))

## POINT Z (1 2 3)

Veamos otros ejemplos:

p <- rbind(c(3.2,4), c(3,4.6), c(3.8,4.4), c(3.5,3.8), c(3.4,3.6), c(3.9,4.5))
(mp <- st_multipoint(p))

## MULTIPOINT ((3.2 4), (3 4.6), (3.8 4.4), (3.5 3.8), (3.4 3.6), (3.9 4.5))

par(mar = c(0.1, 0.1, 1.3, 0.1), mfrow = c(2, 3))
plot(mp, col = 'red')
box()
title("MULTIPOINT")

s1 <- rbind(c(0,3),c(0,4),c(1,5),c(2,5))
(ls <- st_linestring(s1))

## LINESTRING (0 3, 0 4, 1 5, 2 5)

plot(ls, col = 'red')
box()
title("LINESTRING")

s2 <- rbind(c(0.2,3), c(0.2,4), c(1,4.8), c(2,4.8))
s3 <- rbind(c(0,4.4), c(0.6,5))
(mls <- st_multilinestring(list(s1,s2,s3)))

## MULTILINESTRING ((0 3, 0 4, 1 5, 2 5), (0.2 3, 0.2 4, 1 4.8, 2 4.8), (0 4.4, 0.6 5))

plot(mls, col = 'red')
box()
title("MULTILINESTRING")

Los objetos creados se muestra a continuación:

p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2))
plot(pol, border = 'red', col = 'grey', xlim = c(0,4))
box()
title("POLYGON")

p3 <- rbind(c(3,0), c(4,0), c(4,1), c(3,1), c(3,0))
p4 <- rbind(c(3.3,0.3), c(3.8,0.3), c(3.8,0.8), c(3.3,0.8), c(3.3,0.3))[5:1,]
p5 <- rbind(c(3,3), c(4,2), c(4,3), c(3,3))
(mpol <- st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5))))

## MULTIPOLYGON (((0 0, 1 0, 3 2, 2 4, 1 4, 0 0), (1 1, 1 2, 2 2, 1 1)), ((3 0, 4 0, 4 1, 3 1, 3 0), (3.3 0.3, 3.3 0.8, 3.8 0.8, 3.8 0.3, 3.3 0.3)), ((3 3, 4 2, 4 3, 3 3)))

plot(mpol, border = 'red', col = 'grey')
box()
title("MULTIPOLYGON")

(gc <- st_geometrycollection(list(mp, mpol, ls)))

## GEOMETRYCOLLECTION (MULTIPOINT ((3.2 4), (3 4.6), (3.8 4.4), (3.5 3.8), (3.4 3.6), (3.9 4.5)), MULTIPOLYGON (((0 0, 1 0, 3 2, 2 4, 1 4, 0 0), (1 1, 1 2, 2 2, 1 1)), ((3 0, 4 0, 4 1, 3 1, 3 0), (3.3 0.3, 3.3 0.8, 3.8 0.8, 3.8 0.3, 3.3 0.3)), ((3 3, 4 2, 4 3, 3 3))), LINESTRING (0 3, 0 4, 1 5, 2 5))

plot(gc, border = 'grey', col = 'grey')
box()
title("GEOMETRYCOLLECTION")

par(mfrow = c(1, 1))

3.6 Lectura y escritura

Como vimos antes, se indica la dirección del archivo en filename:

filename <- system.file("shape/nc.shp", package="sf")
nc <- st_read(filename)

## Reading layer `nc' from data source 
##   `/Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/sf/shape/nc.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27

Podemos omitir la información de salida con quiet=TRUE o usando una función equivalente read_sf:

nc <- st_read(filename,quiet=TRUE)
nc <- read_sf(filename)

Se puede escribir usando:

st_write(nc, "nc.shp")
write_sf(nc, "nc.shp") # con quiet