3. Organización de los objetos espaciales en R
Esta sección y las siguientes son códigos tomados del libro Análisis espacial con R: Usa R como un Sistema de Información Geográfica, esto se realizó con el fin de practicar para el área de estadística espacial.
Instalación del paquete sf
#install.packages('sf')3.1 Datos vectoriales: modelo simple feature
Cargando la libreria sf
library(sf)3.1.1 Cobertura de puntos
sf maneja datos en dos, tres y hasta cuatro dimensiones.
Las coordenadas se dan en forma de vector o tabla y luego se pasan a una cobertura de puntos de la clase sfg, con la función st_point() para un solo punto o si son varios con st_multipoint()
P <- st_point(c(2, 5)) # punto de dos dimensiones (XY)
class(P)FALSE [1] "XY" "POINT" "sfg"Se observa que P es un punto de la clase sfg, ubicado en el sistema de coordenadas XY.
P <- st_point(c(2, 5, 17, 44), 'XYZM') # Cuatro dimensiones
class(P)FALSE [1] "XYZM" "POINT" "sfg"Se observa que P es un punto de la clase sfg, ubicado en el sistema de coordenadas XYZM.
str(P)FALSE 'XYZM' num [1:4] 2 5 17 44Graficación del punto.
plot(P, axes = T)
Creación de varios puntos.
# Multipoint: coordenadas de varios puntos en 2, 3 o 4 dimensiones
# Crea un vector of coordenadas en x
Xs <- c(2, 4, 5)
# Crea un vector of coordenadas en y
Ys <- c(5, 4, 8)
# Pega Xs y Ys para crear una tabla de coordenadas
coords <- cbind(Xs, Ys)
print(coords)FALSE Xs Ys
FALSE [1,] 2 5
FALSE [2,] 4 4
FALSE [3,] 5 8# Crea el objeto Multipoint (MP)
MP <- st_multipoint(coords)
# Visualización en dos dimensiones.
plot(MP, axes = TRUE)
# Tipo de clase del objeto
class(MP)FALSE [1] "XY" "MULTIPOINT" "sfg"# Se muestra el objeto
print(MP)Se observa que en este caso ya es un multipunto en dos dimensiones.
Multipoint en tres dimensiones
# Multipoint en 3 dimensiones
xyz <- cbind(coords,c(17, 22, 31))
print(xyz)FALSE Xs Ys
FALSE [1,] 2 5 17
FALSE [2,] 4 4 22
FALSE [3,] 5 8 31MP <- st_multipoint(xyz)
print(MP)Un objeto multipunto con tres puntos, en tres dimensiones.
Ahora se mostrara como se pueden junta objetos en colecciones (sfc) con la función st_sfc().
# Colecciones de Simple Features (sfc)
# Crea varios sfg
P1 <- st_point(c(2, 5)); P2 <- st_point(c(4, 4)); P3 <- st_point(c(5, 8))
# Junta varios sfg en un sfc (colección de simple features)
geometria1 <- st_sfc(P1,P2,P3)
# st_sfc(P1,P2, crs = 4326) Proy geográfica LatLong datum WGS84Los objetos anteriores pueden ser más sofisticados si se agrega una tabla con atributos.
# Asociando una geometria sfc con una tabla de atributos (data frame)
# Tabla con ID (campo num) e información adicional (tabla de atributos)
num <- c(1,2,3)
nombre <- c("Pozo","Gasolinera","Pozo")
tabpuntos <- data.frame(cbind(num,nombre))
class(tabpuntos)FALSE [1] "data.frame"print(tabpuntos)FALSE num nombre
FALSE 1 1 Pozo
FALSE 2 2 Gasolinera
FALSE 3 3 PozoUniendo la tabla de atributos con el objeto ‘geometria1’
# sf object
SFP <- st_sf(tabpuntos, geometry = geometria1)
class(SFP) # doble clase: simple feature y dataframeFALSE [1] "sf" "data.frame"print(SFP)FALSE Simple feature collection with 3 features and 2 fields
FALSE Geometry type: POINT
FALSE Dimension: XY
FALSE Bounding box: xmin: 2 ymin: 4 xmax: 5 ymax: 8
FALSE CRS: NA
FALSE num nombre geometry
FALSE 1 1 Pozo POINT (2 5)
FALSE 2 2 Gasolinera POINT (4 4)
FALSE 3 3 Pozo POINT (5 8)Gráficas
plot() despliega el mapa en pantalla: se crea un mapa para cada atributo de la tabla. Si se desea desplegar únicamente la geometría del objeto, en este caso la ubicación de los puntos, se usa la función st_geometry(SFP)
plot(SFP, axes = T)
plot(st_geometry(SFP), axes = T)
Se puede pasar de la clase sf a la clase data.frame, pero en la clase sf de igual manera funcioan los filtros con en la clase data.frame
# Se puede extraer la tabla de atributos de un objeto SFC con
as.data.frame(SFP)FALSE num nombre geometry
FALSE 1 1 Pozo POINT (2 5)
FALSE 2 2 Gasolinera POINT (4 4)
FALSE 3 3 Pozo POINT (5 8)# Selección de elementos dentro de la cobertura
Pozos <- SFP[nombre=="Pozo",]3.1.2. Cobertura de líneas
# Crea 3 objetos "Linestring": simple cadena de coordenadas (vértices)
# Línea L1
X1s <- c(0, 3, 5, 8, 10)
Y1s <- c(0, 3, 4, 8, 10)
Coord1 <- cbind(X1s, Y1s)
# Crea objeto de Clase Linestring
L1 <- st_linestring(Coord1)
class(L1)FALSE [1] "XY" "LINESTRING" "sfg"print(L1)
# Línea L2
X2s <- c(2,1,1)
Y2s <- c(2,4,5)
Coord2 <- cbind(X2s,Y2s)
# Crea objeto de Clase Linestring
L2 <- st_linestring(Coord2)
# Línea 3
X3s <- c(8,8)
Y3s <- c(8,5)
Coord3 <- cbind(X3s,Y3s)
# Crea objeto de Clase Linestring
L3 <- st_linestring(Coord3)
# Crea un objeto Multilineas: conjunto de objetos Linestring
L1L2 <- st_multilinestring(list(L1,L2))
# Junta varios sfg en un sfc (colección de simple features)
geometria2 <- st_sfc(L1,L2,L3)De igual forma como si fueran puntos podemos asociar una tabla de atributos al conjunto de lineas.
# Tabla de atributos
num <- c(1, 2, 3)
code <- c("t", "t", "p")
tipo <- c("Terraceria", "Terraceria", "Pavimentada")
tablineas <- data.frame(cbind(num, tipo, code))
print(tablineas)FALSE num tipo code
FALSE 1 1 Terraceria t
FALSE 2 2 Terraceria t
FALSE 3 3 Pavimentada p# sf object
SFL <- st_sf(tablineas, geometry = geometria2)
plot(SFL, axes=TRUE)
Miremos el conjunto de datos y el tipo de clase
class(SFL)FALSE [1] "sf" "data.frame"print(SFL)FALSE Simple feature collection with 3 features and 3 fields
FALSE Geometry type: LINESTRING
FALSE Dimension: XY
FALSE Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
FALSE CRS: NA
FALSE num tipo code geometry
FALSE 1 1 Terraceria t LINESTRING (0 0, 3 3, 5 4, ...
FALSE 2 2 Terraceria t LINESTRING (2 2, 1 4, 1 5)
FALSE 3 3 Pavimentada p LINESTRING (8 8, 8 5)Al igual que los puntos se pueden filtrar lineas con ciertas características.
# Se puede extraer la tabla de atributos de un SFC con
as.data.frame(SFL)FALSE num tipo code geometry
FALSE 1 1 Terraceria t LINESTRING (0 0, 3 3, 5 4, ...
FALSE 2 2 Terraceria t LINESTRING (2 2, 1 4, 1 5)
FALSE 3 3 Pavimentada p LINESTRING (8 8, 8 5)# Se puede selecionar ciertos rasgos usando la tabla de atributos
print(SFL[tipo == "Pavimentada", ])FALSE Simple feature collection with 1 feature and 3 fields
FALSE Geometry type: LINESTRING
FALSE Dimension: XY
FALSE Bounding box: xmin: 8 ymin: 5 xmax: 8 ymax: 8
FALSE CRS: NA
FALSE num tipo code geometry
FALSE 3 3 Pavimentada p LINESTRING (8 8, 8 5)plot(st_geometry(SFL[tipo=="Terraceria", ]), col = "red", axes = TRUE)
plot(st_geometry(SFL[tipo=="Pavimentada", ]), add = TRUE)
3.1.3. Cobertura de polígonos
Los poligonos son formas cerradas y se crean con la función st_polygon()
La secuencia de los vértices van en el sentido de las manecillas del reloj para delimitar huecos dentro del polígono
La secuencia de los vértices van en el sentido contrario al de las manecillas del reloj para delimitar el contorno externo.
## P1 Polígono forestal al SudEste
## Polygon
# Crea una cadena de coordenadas en X
X1 <- c(5, 10, 10, 6, 5)
# Crea una cadena de coordenadas en Y
# Ojo tiene que cerrar (último par de coord = primero)
Y1 <- c(0, 0, 5, 5, 0)
# Pega X y Y para crear una tabla de coordenadas
c1 <- cbind(X1, Y1)
print(c1)FALSE X1 Y1
FALSE [1,] 5 0
FALSE [2,] 10 0
FALSE [3,] 10 5
FALSE [4,] 6 5
FALSE [5,] 5 0class(c1)FALSE [1] "matrix" "array"Ahora vamos a crear el polígono
# Crea el objeto Polygon. Un polygon es una forma simple cerrada
# eventualmente con hueco(s)
P1 <- st_polygon(list(c1))
plot(P1,axes=T)
Se crearan varios poligonos.
# P2 Polígono forestal al NorOeste
# Crea una cadena de coordenadas en X
X2 <- c(0, 2, 6, 0, 0)
# Crea una cadena de coordenadas en Y
Y2 <- c(4, 4, 10, 10, 4)
# Pega X y Y para crear una tabla de coordenadas
c2 <- cbind(X2,Y2)
# Polígono hueco %%%%%%%%%%%%%% orden coordenadas !!!!!!!!!!!!!
# Crea una cadena de coordenadas en X
X3 <- c(1, 1, 2, 2, 1)
# Crea una cadena de coordenadas en Y
# La 1a y última coordenadas se repiten para "cerrar" el polígono
Y3 <- c(5, 6, 6, 5, 5)
# Pega X y Y para crear una tabla de coordenadas
c3 <- cbind(X3, Y3)
P2 <- st_polygon(list(c2, c3))
plot(P2, axes=T)
Ahora dibujaremos un polígono tipo parcela de agricultura.
# P4 Polígono de agricultura
c3i <- c3[nrow(c3):1, ] # Invierte el orden de las coordenadas
P4 = st_polygon(list(c3i)) # Esta vez no es hueco
# P5 Polígono de área urbana
X5 <- c(0, 5, 6, 10, 10, 6, 2, 0, 0)
Y5 <- c(0, 0, 5, 5, 10, 10, 4, 4, 0)
c5 <- cbind(X5, Y5)
P5 <- st_polygon(list(c5))
plot(P5, axes=T)
Se juntan todos los polígonos anteriores.
# Junta varios sfg en un sfc (colección de simple features)
geometria3 <- st_sfc(P1, P2, P4, P5)Se asocia esta cobertura de polígonos a una tabla de atributos para crear un objeto de la clase sf de forma similar a aquella utilizada para las coberturas de puntos y líneas.
# Tabla de atributos
ID <- c(1, 2, 3, 4)
code <- c("F", "F", "U", "A")
tipo <- c("Bosque", "Bosque", "Urbano", "Agricultura")
tabpol <- data.frame(cbind(ID, code, tipo))
# sf object
SFPol <- st_sf(tabpol, geometry = geometria3)
plot(SFPol, axes=TRUE)
Mirando la clases del dataframe construido
class(SFPol)FALSE [1] "sf" "data.frame"print(SFPol)FALSE Simple feature collection with 4 features and 3 fields
FALSE Geometry type: POLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
FALSE CRS: NA
FALSE ID code tipo geometry
FALSE 1 1 F Bosque POLYGON ((5 0, 10 0, 10 5, ...
FALSE 2 2 F Bosque POLYGON ((0 4, 2 4, 6 10, 0...
FALSE 3 3 U Urbano POLYGON ((1 5, 2 5, 2 6, 1 ...
FALSE 4 4 A Agricultura POLYGON ((0 0, 5 0, 6 5, 10...summary(SFPol)FALSE ID code tipo geometry
FALSE Length:4 Length:4 Length:4 POLYGON:4
FALSE Class :character Class :character Class :character epsg:NA:0
FALSE Mode :character Mode :character Mode :characterLa tabla de atributos permite extraer ciertos rasgos. Las últimas líneas del código a continuación muestran que una colección puede eventualmente juntar diferentes tipos de geometría como puntos, lineas y polígonos.
# Se puede extraer la tabla de atributos de un SFC con
as.data.frame(SFPol)FALSE ID code tipo geometry
FALSE 1 1 F Bosque POLYGON ((5 0, 10 0, 10 5, ...
FALSE 2 2 F Bosque POLYGON ((0 4, 2 4, 6 10, 0...
FALSE 3 3 U Urbano POLYGON ((1 5, 2 5, 2 6, 1 ...
FALSE 4 4 A Agricultura POLYGON ((0 0, 5 0, 6 5, 10...# Se puede selecionar ciertos rasgos usando la tabla de atributos
print(SFPol[tipo=="Bosque", ])FALSE Simple feature collection with 2 features and 3 fields
FALSE Geometry type: POLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
FALSE CRS: NA
FALSE ID code tipo geometry
FALSE 1 1 F Bosque POLYGON ((5 0, 10 0, 10 5, ...
FALSE 2 2 F Bosque POLYGON ((0 4, 2 4, 6 10, 0...Bosque <- SFPol[tipo=="Bosque", ]
plot(st_geometry(SFPol), axes=FALSE)
plot(Bosque, add = TRUE, col = "green")
### Una colección puede eventualmente juntar diferentes tipos de geometría
Detodo <- st_sfc(P1, P2, L1, L2, P1)
plot(Detodo, border = "red", axes = T, col = 'black')
3.2 Datos raster
library(raster)
# Creamos una matriz
m <- matrix(c(1, 2, 3, 4, 2, NA, 2, 2, 3, 3, 3, 1), ncol = 4, nrow = 3, byrow = TRUE)
mFALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 2 3 4
FALSE [2,] 2 NA 2 2
FALSE [3,] 3 3 3 1r <- raster(m)
extent(r) <- extent(c(0, 4, 0, 3))
class(r)FALSE [1] "RasterLayer"
FALSE attr(,"package")
FALSE [1] "raster"print(r)FALSE class : RasterLayer
FALSE dimensions : 3, 4, 12 (nrow, ncol, ncell)
FALSE resolution : 1, 1 (x, y)
FALSE extent : 0, 4, 0, 3 (xmin, xmax, ymin, ymax)
FALSE crs : NA
FALSE source : memory
FALSE names : layer
FALSE values : 1, 4 (min, max)plot(r)
4. Importación/exportación de datos espaciales.
4.1 Importación de archivos shape
shapefile está compuesto por lo menos por tres ficheros:
.shp: almacena las entidades geométricas de los objetos.
.shx: almacena el índice de las entidades geométricas.
.dbf: almacena la tabla de atributos de los objetos.
Es común encontrar un cuarto archivo con extensión .prj que contiene la información referida al sistema de coordenadas.
# Carga el paquete sf
library(sf)
mx <- st_read("recursos-mx/Entidades_latlong.shp")FALSE Reading layer `Entidades_latlong' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\Entidades_latlong.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 32 features and 2 fields
FALSE Geometry type: MULTIPOLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: -118.4076 ymin: 14.5321 xmax: -86.71041 ymax: 32.71865
FALSE Geodetic CRS: WGS 84# Pregunta la clase del objeto espacial
class(mx) # Es un simple feature "sf"FALSE [1] "sf" "data.frame"summary(mx)FALSE NOM_ENT CveEdo geometry
FALSE Length:32 Min. : 1.00 MULTIPOLYGON :32
FALSE Class :character 1st Qu.: 8.75 epsg:4326 : 0
FALSE Mode :character Median :16.50 +proj=long...: 0
FALSE Mean :16.50
FALSE 3rd Qu.:24.25
FALSE Max. :32.00# plotea el mapa (um mapa para cada atributo)
plot(mx, axes = T, cex.axis = 0.8) #cex.axis controla tamaño coordenadas
Si queremos plotear únicamente la geometría
plot(st_geometry(mx))
Conocer el sistema de proyección
st_geometry(mx)FALSE Geometry set for 32 features
FALSE Geometry type: MULTIPOLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: -118.4076 ymin: 14.5321 xmax: -86.71041 ymax: 32.71865
FALSE Geodetic CRS: WGS 84
FALSE First 5 geometries:st_crs(mx)FALSE Coordinate Reference System:
FALSE User input: WGS 84
FALSE wkt:
FALSE GEOGCRS["WGS 84",
FALSE DATUM["World Geodetic System 1984",
FALSE ELLIPSOID["WGS 84",6378137,298.257223563,
FALSE LENGTHUNIT["metre",1]]],
FALSE PRIMEM["Greenwich",0,
FALSE ANGLEUNIT["degree",0.0174532925199433]],
FALSE CS[ellipsoidal,2],
FALSE AXIS["latitude",north,
FALSE ORDER[1],
FALSE ANGLEUNIT["degree",0.0174532925199433]],
FALSE AXIS["longitude",east,
FALSE ORDER[2],
FALSE ANGLEUNIT["degree",0.0174532925199433]],
FALSE ID["EPSG",4326]]Ahora se van a detectar errores de geometría como intersecciones con la función st_is_valid()
carr <- st_read("recursos-mx/carretera.shp")FALSE Reading layer `carretera' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\carretera.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 12424 features and 4 fields
FALSE Geometry type: LINESTRING
FALSE Dimension: XY
FALSE Bounding box: xmin: 1071375 ymin: 324561 xmax: 4076832 ymax: 2349296
FALSE Projected CRS: Lambert_Conformal_Conicst_crs(carr)FALSE Coordinate Reference System:
FALSE User input: Lambert_Conformal_Conic
FALSE wkt:
FALSE PROJCRS["Lambert_Conformal_Conic",
FALSE BASEGEOGCRS["GCS_GRS 1980(IUGG, 1980)",
FALSE DATUM["D_unknown",
FALSE ELLIPSOID["GRS80",6378137,298.257222101,
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]]],
FALSE PRIMEM["Greenwich",0,
FALSE ANGLEUNIT["Degree",0.0174532925199433]]],
FALSE CONVERSION["unnamed",
FALSE METHOD["Lambert Conic Conformal (2SP)",
FALSE ID["EPSG",9802]],
FALSE PARAMETER["Latitude of false origin",12,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8821]],
FALSE PARAMETER["Longitude of false origin",-102,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8822]],
FALSE PARAMETER["Latitude of 1st standard parallel",17.5,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8823]],
FALSE PARAMETER["Latitude of 2nd standard parallel",29.5,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8824]],
FALSE PARAMETER["Easting at false origin",2500000,
FALSE LENGTHUNIT["metre",1],
FALSE ID["EPSG",8826]],
FALSE PARAMETER["Northing at false origin",0,
FALSE LENGTHUNIT["metre",1],
FALSE ID["EPSG",8827]]],
FALSE CS[Cartesian,2],
FALSE AXIS["(E)",east,
FALSE ORDER[1],
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]],
FALSE AXIS["(N)",north,
FALSE ORDER[2],
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]]]No hay información sobre el sistema de coordenadas.
# se define:
st_crs(carr) <- "+proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102
+x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs"
# Test de la validez de los mapas
st_is_valid(mx)FALSE [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE [31] TRUE TRUE4.2 Importación de archivos vector de otros formatos
La función st_read(), que utiliza GDAL, es capaz de importar una gran cantidad de formato vector incluyendo ESRI File Geodatabase (OpenFileGDB) entre otros. Para ver estos formatos, se puede utilizar sf_drivers que permite obtener una lista de los drivers soportados:
drivers_soportados <- st_drivers()
names(drivers_soportados)FALSE [1] "name" "long_name" "write" "copy" "is_raster" "is_vector"
FALSE [7] "vsi"# Lista de los 10 primeros drivers de la lista (hay más!)
head(drivers_soportados[,-2], n = 10)FALSE name write copy is_raster is_vector vsi
FALSE ESRIC ESRIC FALSE FALSE TRUE TRUE TRUE
FALSE FITS FITS TRUE FALSE TRUE TRUE FALSE
FALSE PCIDSK PCIDSK TRUE FALSE TRUE TRUE TRUE
FALSE netCDF netCDF TRUE TRUE TRUE TRUE FALSE
FALSE PDS4 PDS4 TRUE TRUE TRUE TRUE TRUE
FALSE VICAR VICAR TRUE TRUE TRUE TRUE TRUE
FALSE JP2OpenJPEG JP2OpenJPEG FALSE TRUE TRUE TRUE TRUE
FALSE JPEG2000 JPEG2000 FALSE TRUE TRUE TRUE TRUE
FALSE PDF PDF TRUE TRUE TRUE TRUE FALSE
FALSE MBTiles MBTiles TRUE TRUE TRUE TRUE TRUE4.3 Exportación a shape o a otros formatos
Podemos salvar un objeto espacial de la clase sf en formato shape o en otro formato con la función st_write().
# Salva el objeto en formato shape
st_write(mx, "mexico.shp", delete_layer = TRUE)FALSE Deleting layer `mexico' using driver `ESRI Shapefile'
FALSE Writing layer `mexico' to data source `mexico.shp' using driver `ESRI Shapefile'
FALSE Writing 32 features with 2 fields and geometry type Multi Polygon.st_write(mx, dsn = "mx.gpkg", delete_layer = TRUE)FALSE Deleting layer `mx' using driver `GPKG'
FALSE Writing layer `mx' to data source `mx.gpkg' using driver `GPKG'
FALSE Writing 32 features with 2 fields and geometry type Multi Polygon.Dentro de R, es también posible transformar objetos espaciales entre formatos de varios paquetes. Esta opción puede ser particularmente útil para utilizar ciertos paquetes que solo manejan objetos de sp. En sp la función as() permite importar un objeto sf en clase Spatial de sp. Al inverso, st_as_sf() permite pasar de sp a sf.
mx_sp <- as(mx, Class = "Spatial")
# Regreso a sf con st_as_sf():
mx_sf <- st_as_sf(mx_sp, "sf")4.4 Importación / exportación de datos raster
En el paquete raster, la función raster() permite importar imágenes de muchísimos formatos de imagen en objetos RasterLayer. Al contrario, la función writeRaster() permite salvar los objetos RasterLayer en varios formatos de imagen incluyendo ENVI, ESRI, ERDAS, GeoTiff, IDRISI y SAGA (opción format). La opción datatype permite escoger la codificación de los datos. Por ejemplo INT1U, es la codificación en 8 bits con valores de 0 a 255.
library(raster)
# Importa imagen
dem <- raster("recursos-mx/DEM_GTOPO1KM.tif")
class(dem)FALSE [1] "RasterLayer"
FALSE attr(,"package")
FALSE [1] "raster"plot(dem)
extent(dem)FALSE class : Extent
FALSE xmin : -1999500
FALSE xmax : 3000500
FALSE ymin : -3999500
FALSE ymax : 500# Formatos disponibles para salvar
writeFormats()FALSE name long_name
FALSE [1,] "raster" "R-raster"
FALSE [2,] "SAGA" "SAGA GIS"
FALSE [3,] "IDRISI" "IDRISI"
FALSE [4,] "IDRISIold" "IDRISI (img/doc)"
FALSE [5,] "BIL" "Band by Line"
FALSE [6,] "BSQ" "Band Sequential"
FALSE [7,] "BIP" "Band by Pixel"
FALSE [8,] "ascii" "Arc ASCII"
FALSE [9,] "CDF" "NetCDF"
FALSE [10,] "ADRG" "ARC Digitized Raster Graphics"
FALSE [11,] "BAG" "Bathymetry Attributed Grid"
FALSE [12,] "BMP" "MS Windows Device Independent Bitmap"
FALSE [13,] "BT" "VTP .bt (Binary Terrain) 1.3 Format"
FALSE [14,] "BYN" "Natural Resources Canada's Geoid"
FALSE [15,] "CTable2" "CTable2 Datum Grid Shift"
FALSE [16,] "EHdr" "ESRI .hdr Labelled"
FALSE [17,] "ELAS" "ELAS"
FALSE [18,] "ENVI" "ENVI .hdr Labelled"
FALSE [19,] "ERS" "ERMapper .ers Labelled"
FALSE [20,] "FITS" "Flexible Image Transport System"
FALSE [21,] "GPKG" "GeoPackage"
FALSE [22,] "GS7BG" "Golden Software 7 Binary Grid (.grd)"
FALSE [23,] "GSBG" "Golden Software Binary Grid (.grd)"
FALSE [24,] "GTiff" "GeoTIFF"
FALSE [25,] "GTX" "NOAA Vertical Datum .GTX"
FALSE [26,] "HDF4Image" "HDF4 Dataset"
FALSE [27,] "HFA" "Erdas Imagine Images (.img)"
FALSE [28,] "IDA" "Image Data and Analysis"
FALSE [29,] "ILWIS" "ILWIS Raster Map"
FALSE [30,] "INGR" "Intergraph Raster"
FALSE [31,] "ISCE" "ISCE raster"
FALSE [32,] "ISIS2" "USGS Astrogeology ISIS cube (Version 2)"
FALSE [33,] "ISIS3" "USGS Astrogeology ISIS cube (Version 3)"
FALSE [34,] "KRO" "KOLOR Raw"
FALSE [35,] "LAN" "Erdas .LAN/.GIS"
FALSE [36,] "Leveller" "Leveller heightfield"
FALSE [37,] "MBTiles" "MBTiles"
FALSE [38,] "MRF" "Meta Raster Format"
FALSE [39,] "netCDF" "Network Common Data Format"
FALSE [40,] "NGW" "NextGIS Web"
FALSE [41,] "NITF" "National Imagery Transmission Format"
FALSE [42,] "NTv2" "NTv2 Datum Grid Shift"
FALSE [43,] "NWT_GRD" "Northwood Numeric Grid Format .grd/.tab"
FALSE [44,] "PAux" "PCI .aux Labelled"
FALSE [45,] "PCIDSK" "PCIDSK Database File"
FALSE [46,] "PCRaster" "PCRaster Raster File"
FALSE [47,] "PDF" "Geospatial PDF"
FALSE [48,] "PDS4" "NASA Planetary Data System 4"
FALSE [49,] "PNM" "Portable Pixmap Format (netpbm)"
FALSE [50,] "RMF" "Raster Matrix Format"
FALSE [51,] "ROI_PAC" "ROI_PAC raster"
FALSE [52,] "RRASTER" "R Raster"
FALSE [53,] "RST" "Idrisi Raster A.1"
FALSE [54,] "SAGA" "SAGA GIS Binary Grid (.sdat, .sg-grd-z)"
FALSE [55,] "SGI" "SGI Image File Format 1.0"
FALSE [56,] "Terragen" "Terragen heightfield"
FALSE [57,] "VICAR" "MIPL VICAR file"Vamos a guardar el raster en otro tipo de formato
# Salva el raster en formato LAN
writeRaster(dem, filename = "DEM_Mx.lan", format = "LAN", overwrite = TRUE, datatype = "INT2S")5 Operaciones básicas de SIG (vector)
5.1 Algunas operaciones de análisis espacial
Algunas funciones de sf permiten verificar la geometría de objetos y calcular la longitud de líneas simples y el área de polígonos. La segunda parte del código, nos permite crear un mapa con los objetos espaciales que vamos a analizar. Los números indican el número de los puntos y las líneas.
library(sf)
# Importación de los mapitas con st_read (paquete sf)
# Mapa puntos
SFP <- st_read("recursos-mx/SFP.shp")FALSE Reading layer `SFP' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\SFP.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 3 features and 2 fields
FALSE Geometry type: POINT
FALSE Dimension: XY
FALSE Bounding box: xmin: 2 ymin: 4 xmax: 5 ymax: 8
FALSE CRS: NA# Mapa de líneas
SFL <- st_read("recursos-mx/SFL.shp")FALSE Reading layer `SFL' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\SFL.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 3 features and 3 fields
FALSE Geometry type: LINESTRING
FALSE Dimension: XY
FALSE Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
FALSE CRS: NA# Mapa de polígonos
SFPol <- st_read("recursos-mx/SFPol.shp")FALSE Reading layer `SFPol' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\SFPol.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 4 features and 3 fields
FALSE Geometry type: POLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
FALSE CRS: NA# Test de la validez de polígonos
st_is_valid(SFPol)FALSE [1] TRUE TRUE TRUE TRUE# Test si líneas son simples
st_is_simple(SFL)FALSE [1] TRUE TRUE TRUE# Calcula área de polígonos
st_area(SFPol)FALSE [1] 22.5 23.0 1.0 53.5# Calcula longitud líneas (no funcion para multiline)
st_length(SFL)FALSE [1] 14.307136 3.236068 3.000000# Mapita de los datos
plot(st_geometry(SFPol["ID"]), col = c("Dark Green", "green", "grey", "white"))
plot(st_geometry(SFP), add = TRUE)
plot(st_geometry(SFP) + c(0.5, 0), add = TRUE, pch = c(49, 50, 51), cex = 1)
plot(st_geometry(SFL), lty = 2, lwd = 2, col = "red", add = T)
# agrega el número de las líneas
text(5, 5.5, "1", pos = 4, col = "red", cex = 1)
text(1.2, 3, "2", pos = 4, col = "red", cex = 1)
text(7.8, 6, "3", pos = 4, col = "red", cex = 1)
En un primer paso, vamos a utilizar la función st_intersects(), que permite probar si los elementos de objetos espaciales se intersectan o no. La opción sparse permite manejar dos formas de presentación de los resultados: como una matriz (sparse = FALSE), en este case de 3 x 4 (3 líneas, 4 polígonos) o de una forma condensada (sparse = TRUE). st_distance() permite calcular la distancia euclidiana más corta entre elementos de dos objetos y presenta también los resultados en forma de matriz.
## Operadores lógicos binarios
# Intersección de puntos y polígonos
st_intersects(SFP, SFPol, sparse = FALSE) # sparse = FALSE Forma condesadaFALSE [,1] [,2] [,3] [,4]
FALSE [1,] FALSE TRUE TRUE FALSE
FALSE [2,] FALSE FALSE FALSE TRUE
FALSE [3,] FALSE FALSE FALSE TRUEst_intersects(SFP, SFPol, sparse = TRUE) # spase = TRUE Forma matricialFALSE Sparse geometry binary predicate list of length 3, where the predicate
FALSE was `intersects'
FALSE 1: 2, 3
FALSE 2: 4
FALSE 3: 4# Intersección de líneas y polígonos
st_intersects(SFL, SFPol, sparse = FALSE)FALSE [,1] [,2] [,3] [,4]
FALSE [1,] FALSE FALSE FALSE TRUE
FALSE [2,] FALSE TRUE TRUE TRUE
FALSE [3,] TRUE FALSE FALSE TRUEst_intersects(SFL, SFPol, sparse = TRUE)FALSE Sparse geometry binary predicate list of length 3, where the predicate
FALSE was `intersects'
FALSE 1: 4
FALSE 2: 2, 3, 4
FALSE 3: 1, 4# Distancia más corta entre elementos de dos objetos
st_distance(SFP, SFPol)FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 3.922323 0.0000000 0.000000 0.5547002
FALSE [2,] 1.765045 1.6641006 2.236068 0.0000000
FALSE [3,] 3.162278 0.2773501 3.605551 0.0000000Finalmente, st_intersection() permite realizar una intersección geométrica entre dos objetos, obteniendo como resultado un objeto espacial cuyas geometría y tabla de atributos recoge la información de ambos objetos donde estos se sobrelapan. En el código a continuación, esta función permite calcular la longitud de los segmentos que coinciden espacialmente con las diferentes categorías de los polígonos.
# Operadores geométricos (intersección geométrica)
SFPxSFPol <- st_intersection(SFP,SFPol)
print(SFPxSFPol) # punto 1 es duplicado por pertenecer a 2 polígonosFALSE Simple feature collection with 4 features and 5 fields
FALSE Geometry type: POINT
FALSE Dimension: XY
FALSE Bounding box: xmin: 2 ymin: 4 xmax: 5 ymax: 8
FALSE CRS: NA
FALSE num nombre ID code tipo geometry
FALSE 1 1 Pozo 2 F Bosque POINT (2 5)
FALSE 1.1 1 Pozo 3 U Urbano POINT (2 5)
FALSE 2 2 Gasolinera 4 A Agricultura POINT (4 4)
FALSE 3 3 Pozo 4 A Agricultura POINT (5 8)SFLxSFPol <- st_intersection(SFL,SFPol)
print(SFLxSFPol) # hay puntos y líneasFALSE Simple feature collection with 6 features and 6 fields
FALSE Geometry type: GEOMETRY
FALSE Dimension: XY
FALSE Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
FALSE CRS: NA
FALSE num tipo code ID code.1 tipo.1 geometry
FALSE 3 3 Pavimentada p 1 F Bosque POINT (8 5)
FALSE 2 2 Terraceria t 2 F Bosque LINESTRING (1 4, 1 5)
FALSE 2.1 2 Terraceria t 3 U Urbano POINT (1 5)
FALSE 1 1 Terraceria t 4 A Agricultura LINESTRING (0 0, 3 3, 5 4, ...
FALSE 2.2 2 Terraceria t 4 A Agricultura LINESTRING (2 2, 1 4)
FALSE 3.1 3 Pavimentada p 4 A Agricultura LINESTRING (8 8, 8 5)plot(SFLxSFPol["tipo.1"],lwd=3)
# extraemos solo las líneas
SFLxSFPol_l <- st_collection_extract(SFLxSFPol,type="LINESTRING")
longitud <- st_length(SFLxSFPol_l) # Cálculo longitud
SFLxSFPol_l$long <- longitud # Resultados en nueva columna de la tabla
# Suma la longitud de los segmentos de cada tipo de cubierta
Suma_l <- aggregate(long ~ tipo.1, FUN = sum, data = SFLxSFPol_l)
print(Suma_l)FALSE tipo.1 long
FALSE 1 Agricultura 19.5432
FALSE 2 Bosque 1.0000Ahora haciendolo de una manera más condesada con el operador pipe
# Con pipe
longitud <- st_intersection(SFL,SFPol) %>% st_collection_extract(type = "LINESTRING") %>% st_length()
print(longitud)FALSE [1] 1.000000 14.307136 2.236068 3.0000005.2 Análisis espacial en formato vector
En este sección, vamos a realizar un análisis basado en la intersección entre dos mapas. Vamos a importar dos mapas en formato shape: Un mapa de los Estados Mexicanos y otro de las principales carreteras. Vamos a visualizar a tabla de atributos de estos mapas y reproyectar el mapa de los Estados que está en Lat/Long a Cónica Conforme de Lambert que es la proyección en la cual se encuentra el mapa de carreteras. Finalmente, calculamos el área de cada estado. Es importante notar que dividimos el área obtenida (en \(m^2\)) por un millón para obtener \(km^22\) pero que el sistema indica \(m^2\) como unidades.
# Carga los paquetes
library(sf)
library(raster)
# Importación de los mapas con st_read (paquete sf)
# Mapa de los estados de México
mx <- st_read("recursos-mx/Entidades_latlong.shp")FALSE Reading layer `Entidades_latlong' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\Entidades_latlong.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 32 features and 2 fields
FALSE Geometry type: MULTIPOLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: -118.4076 ymin: 14.5321 xmax: -86.71041 ymax: 32.71865
FALSE Geodetic CRS: WGS 84st_crs(mx) # sistema de coordenadasFALSE Coordinate Reference System:
FALSE User input: WGS 84
FALSE wkt:
FALSE GEOGCRS["WGS 84",
FALSE DATUM["World Geodetic System 1984",
FALSE ELLIPSOID["WGS 84",6378137,298.257223563,
FALSE LENGTHUNIT["metre",1]]],
FALSE PRIMEM["Greenwich",0,
FALSE ANGLEUNIT["degree",0.0174532925199433]],
FALSE CS[ellipsoidal,2],
FALSE AXIS["latitude",north,
FALSE ORDER[1],
FALSE ANGLEUNIT["degree",0.0174532925199433]],
FALSE AXIS["longitude",east,
FALSE ORDER[2],
FALSE ANGLEUNIT["degree",0.0174532925199433]],
FALSE ID["EPSG",4326]]summary(mx) # resumen descriptivoFALSE NOM_ENT CveEdo geometry
FALSE Length:32 Min. : 1.00 MULTIPOLYGON :32
FALSE Class :character 1st Qu.: 8.75 epsg:4326 : 0
FALSE Mode :character Median :16.50 +proj=long...: 0
FALSE Mean :16.50
FALSE 3rd Qu.:24.25
FALSE Max. :32.00# Mapa de carreteras
carr <- st_read("recursos-mx/carretera.shp")FALSE Reading layer `carretera' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\carretera.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 12424 features and 4 fields
FALSE Geometry type: LINESTRING
FALSE Dimension: XY
FALSE Bounding box: xmin: 1071375 ymin: 324561 xmax: 4076832 ymax: 2349296
FALSE Projected CRS: Lambert_Conformal_Conic#Sistema de coordenadas
st_crs(carr)FALSE Coordinate Reference System:
FALSE User input: Lambert_Conformal_Conic
FALSE wkt:
FALSE PROJCRS["Lambert_Conformal_Conic",
FALSE BASEGEOGCRS["GCS_GRS 1980(IUGG, 1980)",
FALSE DATUM["D_unknown",
FALSE ELLIPSOID["GRS80",6378137,298.257222101,
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]]],
FALSE PRIMEM["Greenwich",0,
FALSE ANGLEUNIT["Degree",0.0174532925199433]]],
FALSE CONVERSION["unnamed",
FALSE METHOD["Lambert Conic Conformal (2SP)",
FALSE ID["EPSG",9802]],
FALSE PARAMETER["Latitude of false origin",12,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8821]],
FALSE PARAMETER["Longitude of false origin",-102,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8822]],
FALSE PARAMETER["Latitude of 1st standard parallel",17.5,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8823]],
FALSE PARAMETER["Latitude of 2nd standard parallel",29.5,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8824]],
FALSE PARAMETER["Easting at false origin",2500000,
FALSE LENGTHUNIT["metre",1],
FALSE ID["EPSG",8826]],
FALSE PARAMETER["Northing at false origin",0,
FALSE LENGTHUNIT["metre",1],
FALSE ID["EPSG",8827]]],
FALSE CS[Cartesian,2],
FALSE AXIS["(E)",east,
FALSE ORDER[1],
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]],
FALSE AXIS["(N)",north,
FALSE ORDER[2],
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]]]# Atributos del mapa de carretera
names(carr)FALSE [1] "ADMINISTRA" "ENTIDAD" "NO_CARRILE" "SHAPE_len" "geometry"# Resumen descriptivo
summary(carr)FALSE ADMINISTRA ENTIDAD NO_CARRILE SHAPE_len
FALSE Length:12424 Length:12424 Length:12424 Min. : 2.23
FALSE Class :character Class :character Class :character 1st Qu.: 1346.58
FALSE Mode :character Mode :character Mode :character Median : 3780.05
FALSE Mean : 6843.31
FALSE 3rd Qu.: 8942.17
FALSE Max. :118873.89
FALSE geometry
FALSE LINESTRING :12424
FALSE epsg:NA : 0
FALSE +proj=lcc ...: 0
FALSE
FALSE
FALSE # Reproyecta mx a LCC
mx_lcc <- st_transform(mx, st_crs(carr))
st_crs(mx_lcc)FALSE Coordinate Reference System:
FALSE User input: Lambert_Conformal_Conic
FALSE wkt:
FALSE PROJCRS["Lambert_Conformal_Conic",
FALSE BASEGEOGCRS["GCS_GRS 1980(IUGG, 1980)",
FALSE DATUM["D_unknown",
FALSE ELLIPSOID["GRS80",6378137,298.257222101,
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]]],
FALSE PRIMEM["Greenwich",0,
FALSE ANGLEUNIT["Degree",0.0174532925199433]]],
FALSE CONVERSION["unnamed",
FALSE METHOD["Lambert Conic Conformal (2SP)",
FALSE ID["EPSG",9802]],
FALSE PARAMETER["Latitude of false origin",12,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8821]],
FALSE PARAMETER["Longitude of false origin",-102,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8822]],
FALSE PARAMETER["Latitude of 1st standard parallel",17.5,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8823]],
FALSE PARAMETER["Latitude of 2nd standard parallel",29.5,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8824]],
FALSE PARAMETER["Easting at false origin",2500000,
FALSE LENGTHUNIT["metre",1],
FALSE ID["EPSG",8826]],
FALSE PARAMETER["Northing at false origin",0,
FALSE LENGTHUNIT["metre",1],
FALSE ID["EPSG",8827]]],
FALSE CS[Cartesian,2],
FALSE AXIS["(E)",east,
FALSE ORDER[1],
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]],
FALSE AXIS["(N)",north,
FALSE ORDER[2],
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]]]# Calcula el área de los Estados
mx_lcc$AreaEdo <- st_area(mx_lcc) / 1000000 # km2
# Plotea los estados y carreteras juntos
plot(st_geometry(carr), col = "blue", axes = TRUE)
plot(st_geometry(mx_lcc), border = "red", add = TRUE)
A continuación, vamos a calcular la longitud de carreteras federales en cada Estado del sureste de México (Tabasco, Campeche, Quintana Roo y Yucatán). En sf, la selección de ciertos rasgos de un mapa se realiza con base en la tabla de atributos. Por ejemplo, carr[ADMINISTRA=“Federal”,] nos permite selecionar las líneas que corresponden a carreteras federales. Debido a que el mapa de carreteras ”carr” es un objeto sf, el resutado de la selección es también un objeto espacial sf y no una simple tabla. La función st_intersection() nos permite realizar la intersección entre los dos mapas. En este caso, el objeto obtenido son líneas (carreteras federales). La información del Estado en la cual se encuentran los segmentos de carretera está en la tabla de atributos. Se calcula la longitud de cada segmento, y se suman los valores de longitud para cada Estado.
# Qué longitud de carreteras federales en c/ estado del sureste?
# Selección carr federales
carrFed <- carr[ADMINISTRA = "Federal", ]
# Selección estados del sureste (Tabasco, Campeche, Quintana Roo y Yucatán)
sureste <- mx_lcc[mx_lcc$NOM_ENT %in% c("Tabasco", "Campeche",
"Quintana Roo", "Yucatan"),]
class(sureste)FALSE [1] "sf" "data.frame"plot(st_geometry(sureste))
# Intersección entre mapas de carr federales y Edos sureste
carrFedXedos <- st_intersection(carrFed, sureste)
names(carrFedXedos)FALSE [1] "ADMINISTRA" "ENTIDAD" "NO_CARRILE" "SHAPE_len" "NOM_ENT"
FALSE [6] "CveEdo" "AreaEdo" "geometry"summary(carrFedXedos)FALSE ADMINISTRA ENTIDAD NO_CARRILE SHAPE_len
FALSE Length:1177 Length:1177 Length:1177 Min. : 2.23
FALSE Class :character Class :character Class :character 1st Qu.: 1820.08
FALSE Mode :character Mode :character Mode :character Median : 5673.34
FALSE Mean : 8697.05
FALSE 3rd Qu.: 11550.34
FALSE Max. :118873.89
FALSE NOM_ENT CveEdo AreaEdo geometry
FALSE Length:1177 Min. : 4.00 Min. :24695 LINESTRING :1147
FALSE Class :character 1st Qu.:23.00 1st Qu.:24695 MULTILINESTRING: 30
FALSE Mode :character Median :27.00 Median :39533 epsg:NA : 0
FALSE Mean :24.38 Mean :38158 +proj=lcc ... : 0
FALSE 3rd Qu.:31.00 3rd Qu.:44556
FALSE Max. :31.00 Max. :57277# Convierte los objetos "Multilines" en "Lines" (para cálculo longitud)
carrFedXedos_l <- st_collection_extract(carrFedXedos, type = "LINESTRING")
long <- st_length(carrFedXedos_l) / 1000 # km
head(long)FALSE Units: [m]
FALSE [1] 1.219532 1.133280 19.162872 14.348592 1.500849 4.092854# Agrega columna long a la tabla de atributos
carrFedXedos_l$long <- long
# Suma la longitud de los segmentos de cada tipo de cubierta
Suma_long <- aggregate(long ~ NOM_ENT, FUN = sum, data = carrFedXedos_l)
print(Suma_long)FALSE NOM_ENT long
FALSE 1 Campeche 2060.536 [m]
FALSE 2 Quintana Roo 1618.576 [m]
FALSE 3 Tabasco 2382.986 [m]
FALSE 4 Yucatan 3738.168 [m]En seguida, se calcula la proporción del territorio de la región sureste que se encuentra a menos de 20 km de una carretera federal. Para ello, se elabora un área buffer de 20,000 m alrededor de las carreteras. Las etapas siguientes son similares al cálculo realizado para la longitud de carretera. Se puede observar que la mayoría de las operaciones son operaciones tabulares que se hacen de la misma forma que con cualquier tabla dataframe.
# Que proporción del territorio estatal está a menos de 20 km de una
# carretera federal?
# Creación de un buffer de 20 km alrededor carreteras fed
buf <- st_buffer(carrFed, dist = 20000)
plot(st_geometry(buf), axes = F)
# Intersección entre buffer de carr federales y Edos sureste
carrbufXedos <- st_intersection(buf, sureste)
area <- st_area(carrbufXedos) / 1000000 # km2
head(area) # Ojo: Pone m2 por default, en este caso son km2FALSE Units: [m^2]
FALSE [1] 38.39414 274.58121 572.58037 680.35825 1197.92400 1329.92446carrbufXedos$area <- area
# Suma la longitud de los segmentos de cada tipo de cubierta
Suma_area <- aggregate(area ~ NOM_ENT, FUN = sum, data = carrbufXedos)
print(Suma_area)FALSE NOM_ENT area
FALSE 1 Campeche 276058.4 [m^2]
FALSE 2 Quintana Roo 176800.9 [m^2]
FALSE 3 Tabasco 506812.3 [m^2]
FALSE 4 Yucatan 704789.8 [m^2]# junta a sureste la tabla Suma_area (por clave en común)
sureste <- merge(sureste, Suma_area)
# calcula la proporción del área del buffer / área del estado
sureste$prop <- sureste$area / sureste$AreaEdoEn el código a continuación, se relaciona la tabla de abributo del mapa de los Estados con una tabla externa (población de los estados) y se calcula un índice basado en variables representadas por diferente columnas de la tabla. Finalmente, se exporta el objeto espacial obtenido a formato shape
# Población Censo de Población y Vivienda 2010
tab_pop <- read.csv('recursos-mx/Pop2010_INEGI.csv')
head(tab_pop)FALSE NumEdo Estado Poptotal Hombre Mujer
FALSE 1 1 Aguascalientes 1184996 576638 608358
FALSE 2 2 Baja California 3155070 1591610 1563460
FALSE 3 3 Baja California Sur 637026 325433 311593
FALSE 4 4 Campeche 822441 407721 41472
FALSE 5 5 Coahuila de Zaragoza 2748391 1364197 1384194
FALSE 6 6 Colima 650555 32279 327765head(mx_lcc)FALSE Simple feature collection with 6 features and 3 fields
FALSE Geometry type: MULTIPOLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: 907821.8 ymin: 485730.8 xmax: 2925387 ymax: 2349609
FALSE Projected CRS: Lambert_Conformal_Conic
FALSE NOM_ENT CveEdo geometry AreaEdo
FALSE 1 Distrito Federal 9 MULTIPOLYGON (((2802176 843... 1486.481 [m^2]
FALSE 2 Guerrero 12 MULTIPOLYGON (((2723198 539... 63565.113 [m^2]
FALSE 3 Mexico 15 MULTIPOLYGON (((2717219 921... 22226.643 [m^2]
FALSE 4 Morelos 17 MULTIPOLYGON (((2808476 786... 4859.432 [m^2]
FALSE 5 Sinaloa 25 MULTIPOLYGON (((2050677 124... 56802.974 [m^2]
FALSE 6 Baja California 2 MULTIPOLYGON (((1458026 185... 73535.910 [m^2]# Junta las dos tablas con la clave numérica
mx_lcc <- merge(mx_lcc, tab_pop, by.x = "CveEdo", by.y = "NumEdo")
head(mx_lcc)FALSE Simple feature collection with 6 features and 7 fields
FALSE Geometry type: MULTIPOLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: 907821.8 ymin: 692158.7 xmax: 3859532 ymax: 2349609
FALSE Projected CRS: Lambert_Conformal_Conic
FALSE CveEdo NOM_ENT AreaEdo Estado Poptotal
FALSE 1 1 Aguascalientes 5558.679 [m^2] Aguascalientes 1184996
FALSE 2 2 Baja California 73535.910 [m^2] Baja California 3155070
FALSE 3 3 Baja California Sur 73986.177 [m^2] Baja California Sur 637026
FALSE 4 4 Campeche 57277.175 [m^2] Campeche 822441
FALSE 5 5 Coahuila de Zaragoza 150670.960 [m^2] Coahuila de Zaragoza 2748391
FALSE 6 6 Colima 5754.971 [m^2] Colima 650555
FALSE Hombre Mujer geometry
FALSE 1 576638 608358 MULTIPOLYGON (((2470518 115...
FALSE 2 1591610 1563460 MULTIPOLYGON (((1458026 185...
FALSE 3 325433 311593 MULTIPOLYGON (((1694646 122...
FALSE 4 407721 41472 MULTIPOLYGON (((3544897 946...
FALSE 5 1364197 1384194 MULTIPOLYGON (((2469954 197...
FALSE 6 32279 327765 MULTIPOLYGON (((1158659 767...names(mx_lcc)FALSE [1] "CveEdo" "NOM_ENT" "AreaEdo" "Estado" "Poptotal" "Hombre" "Mujer"
FALSE [8] "geometry"# Densidad de población (núm hbs / km2)
mx_lcc$dens <- mx_lcc$Poptotal/mx_lcc$AreaEdo
# El índice de masculinidad, también llamado razón de sexo es un índice
# demográfico que expresa la razón de hombres por mujeres en un
# determinado territorio, expresado # por lo tanto en %. Se calcula
# usando la fórmula: 100 * H/M
mx_lcc$IM <- 100 * mx_lcc$Hombre/mx_lcc$Mujer
# salva el objeto en shape
st_write(mx_lcc,"mx_lcc.shp",delete_layer = TRUE)FALSE Deleting layer `mx_lcc' using driver `ESRI Shapefile'
FALSE Writing layer `mx_lcc' to data source `mx_lcc.shp' using driver `ESRI Shapefile'
FALSE Writing 32 features with 9 fields and geometry type Multi Polygon.6. Operaciones básicas de SIG (raster)
6.1 Algunas operaciones de análisis espacial
Presentamos aquí algunas funciones de análisis espacial del paquete raster, que son más fáciles de entender con datos ”minimalistas”. En la segunda sección del capítulo, llevamos a cabo algunas operaciones con mapas de uso y cubierta del suelo del Estado de Michoacán, México. Para una revisión más profunda del paquete raster, consultar Hijmans (2017).
Vamos a crear algunos datos raster de dimensión muy reducida (3 x 4 celdas) para mostrar fácilmente los resultados obtenidos, observando los mapas raster como matrices. Los datos raster puede eventualmente manejarse como imágenes multibanda con stack() o brick(). Los comandos res(), dim() y xmax() permiten obtener algunas de las características de los mapas como resolución (tamaño de las celdas), dimensión (número de filas, columnas y bandas) y coordenadas extremas. cellStats() permite calcular índices sobre el conjunto de la imagen (tomando en cuenta todas las celdas).
library(raster) # se llama la librería
############## Datos muy sencillos
# Creamos unas matrices
m1 <- matrix(c(1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 3, 3), ncol = 4, nrow = 3, byrow = TRUE)
m2 <- matrix(c(1, 2, 3, 4, 2, NA, 2, 2, 3, 3, 3, 1), ncol = 4, nrow = 3, byrow = TRUE)
print(m1) # mostramos por consola las matrizFALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 1 2 2
FALSE [2,] 1 1 2 2
FALSE [3,] 1 1 3 3print(m2) # mostramos por consola la matrizFALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 2 3 4
FALSE [2,] 2 NA 2 2
FALSE [3,] 3 3 3 1r1 <- raster(m1)
r2 <- raster(m2)
extent(r1) <- extent(r2) <- extent(c(0, 4, 0, 3))
### stack
colortable(r1) <- c("red", "blue", "green")
colortable(r2)<- c("red", "blue", "green", "yellow")
r12 <- stack(r1, r2)
plot(r1, axes = TRUE)
plot(r2, axes = TRUE)
Algunas caracteristicas del raster
# Características del raster
res(r1) # ResoluciónFALSE [1] 1 1dim(r12) # dimensiónFALSE [1] 3 4 2xmax(r1) # existen también xmax, ymin y ymaxFALSE [1] 4cellStats(r1,"sum") # suma de todos los pixelesFALSE [1] 20cellStats(r1,"sd") # desv estándar de todos los pixelesFALSE [1] 0.7784989El paquete raster permite realizar todas las operaciones de un SIG. Las operaciones de álgebra de mapa se aplican pixel por pixel. La función overlay() permite utilizar funciones definidas por el usuario.
# Algebra de mapa
sum1 <- r1 + 2 # suma una constante a cada componente.
print(as.matrix(sum1))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 3 3 4 4
FALSE [2,] 3 3 4 4
FALSE [3,] 3 3 5 5sum2 <- r1 + 2*r2
print(as.matrix(sum2))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 3 5 8 10
FALSE [2,] 5 NA 6 6
FALSE [3,] 7 7 9 5sum3 <- overlay(r1, r2, sum1, fun=function(x, y, z){ x + 2*y -z} )
print(as.matrix(sum3))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 0 2 4 6
FALSE [2,] 2 NA 2 2
FALSE [3,] 4 4 4 0Es posible ”pegar” entre ellos varios mapas para elaborar mosaicos con merge(). Al contrario, crop() permite recortar una imagen con base en coordenadas extremas. Estas operaciones se basan en el extent de los mapas, que contiene las coordenadas extremas (esquinas inferior izquierda y superior derecha) en el orden xmin, xmax, ymin e ymax. La figura 4 ilustra las operaciones de elaboración y recorte de un mosaico realizadas en el código a continuación.
# Operaciones tipo SIG
# Mosaico (merge) y recorte (crop)
# Extent de r1 (xmin, xmin, ymin, ymax)
extent(r1)FALSE class : Extent
FALSE xmin : 0
FALSE xmax : 4
FALSE ymin : 0
FALSE ymax : 3m3 <- matrix(c(5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5), ncol = 4, nrow = 3, byrow = TRUE)
print(m3)FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 5 5 5 5
FALSE [2,] 5 5 5 5
FALSE [3,] 5 5 5 5r3 <- raster(m3)
# extent(xmin, xmax, ymin, ymax)
extent(r3) <- extent(c(0, 4, 3, 6)) # Está al norte de r1
mosaico <- merge(r1, r3)
extent(mosaico)FALSE class : Extent
FALSE xmin : 0
FALSE xmax : 4
FALSE ymin : 0
FALSE ymax : 6print(as.matrix(mosaico))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 5 5 5 5
FALSE [2,] 5 5 5 5
FALSE [3,] 5 5 5 5
FALSE [4,] 1 1 2 2
FALSE [5,] 1 1 2 2
FALSE [6,] 1 1 3 3extent_corte <- extent(c(1, 3, 2, 4))
corte <- crop(mosaico, extent_corte) # corte de raster
print(as.matrix(corte))FALSE [,1] [,2]
FALSE [1,] 5 5
FALSE [2,] 1 2A continuación, presentamos algunas de las maneras de reclasificar los valores de un mapa disponibles en el paquete raster. Con r2[r2 > 2] <- 6, estamos asignando el valor 6 a las celdas de r2 que tienen un valor estrictamente superior a 2. Para reclasificaciones que involucran varios intervalos, es más fácil utilizar tablas de reclasificación. Estas tablas son matrices con tres columnas: las dos primeras indican los umbrales del intervalo de reclasificación y la tercera el valor de salida. El primer valor no está incluido en el intervalo, el segundo si. Por ejemplo, la primera fila de la tabla de reclasificación 0 2 0 indica que los valores estrictamente superior a cero e inferior o igual a 2 se reclasifican en cero. La segunda fila 2 5 1 permite reclasificar los valores superiores a 2 hasta 5 (incluido) en el valor uno.
# Reclasificaciones
print(as.matrix(r2))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 2 3 4
FALSE [2,] 2 NA 2 2
FALSE [3,] 3 3 3 1r2[r2 > 2] <- 6
print(as.matrix(r2))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 2 6 6
FALSE [2,] 2 NA 2 2
FALSE [3,] 6 6 6 1# Tabla de reclasificación (rangos)
m <- c(0, 2, 0, 2, 5, 1)
tabla_reclas <- matrix(m, ncol = 3, byrow = TRUE)
print(tabla_reclas)FALSE [,1] [,2] [,3]
FALSE [1,] 0 2 0
FALSE [2,] 2 5 1reclas <- reclassify(mosaico, tabla_reclas)
print(as.matrix(reclas))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 1 1 1
FALSE [2,] 1 1 1 1
FALSE [3,] 1 1 1 1
FALSE [4,] 0 0 0 0
FALSE [5,] 0 0 0 0
FALSE [6,] 0 0 1 1# Substitución de valores
print(as.matrix(r2))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 2 6 6
FALSE [2,] 2 NA 2 2
FALSE [3,] 6 6 6 1r2[r2 == 6] <- 9
print(as.matrix(r2))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 2 9 9
FALSE [2,] 2 NA 2 2
FALSE [3,] 9 9 9 1# con tabla
tab_subs <- data.frame(id = c(1, 2, 3), v = c(1, 56, 84))
print(tab_subs)FALSE id v
FALSE 1 1 1
FALSE 2 2 56
FALSE 3 3 84sub <- subs(r1, tab_subs)
print(as.matrix(r1))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 1 2 2
FALSE [2,] 1 1 2 2
FALSE [3,] 1 1 3 3print(as.matrix(sub))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 1 56 56
FALSE [2,] 1 1 56 56
FALSE [3,] 1 1 84 84Se puede modificar el arreglo espacial de las celdas agrupando celdas. La función aggregate() permite reagrupar varios pixeles, obteniendo un mapa raster de menor resolución espacial. La función fun permite controlar el cálculo del valor del pixel de baja resolución que corresponde a varias celdas en los datos de entrada. En los dos ejemplos a continuación se reagrupan 2x2 pixeles. En la primera opción se suman los valores de las cuatro celdas, en la segunda se escoge el valor más común (mayoritario).
# Remuestreos
# Agregación espacial
agreg <- aggregate(mosaico, fact = 2, fun = "sum")
print(as.matrix(mosaico))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 5 5 5 5
FALSE [2,] 5 5 5 5
FALSE [3,] 5 5 5 5
FALSE [4,] 1 1 2 2
FALSE [5,] 1 1 2 2
FALSE [6,] 1 1 3 3print(as.matrix(agreg))FALSE [,1] [,2]
FALSE [1,] 20 20
FALSE [2,] 12 14
FALSE [3,] 4 10agreg2 <- aggregate(mosaico, fact = 2, fun = modal, na.rm = TRUE)
print(as.matrix(agreg2))FALSE [,1] [,2]
FALSE [1,] 5 5
FALSE [2,] 5 5
FALSE [3,] 1 2El paquete raster permite también realizar remuestreo con los métodos comúnmente utilizados en procesamiento de imágenes y en percepción remota.
# Remuestreo
# Imagen de mismo tamaño que mosaico con menos filas
m4 <- matrix(rep(1, 20), ncol = 4, nrow = 5)
r4 <- raster(m4)
extent(r4) <- extent(mosaico)
resv <- resample(mosaico, r4, method = "ngb") # vecino más cercano
resb <- resample(mosaico, r4, method = "bilinear") # bilinear
print(as.matrix(resv))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 5 5 5 5
FALSE [2,] 5 5 5 5
FALSE [3,] 1 1 2 2
FALSE [4,] 1 1 2 2
FALSE [5,] 1 1 3 3print(as.matrix(resb))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 5 5 5.0 5.0
FALSE [2,] 5 5 5.0 5.0
FALSE [3,] 3 3 3.5 3.5
FALSE [4,] 1 1 2.0 2.0
FALSE [5,] 1 1 2.9 2.9El paquete raster permite realizar operaciones de filtrado espacial conocido como convolución (kernel). Las dos operaciones focales a continuación son equivalentes y calculan el promedio del valores de las celdas en una ventana móvil de 3 x 3 pixeles. El cálculo se realiza únicamente para las celdas que tienen ocho vecinos. Sin embargo, la opción ”pad = TRUE” permite evitar este efecto en los bordes de la imagen (ver la ayuda).
# Operaciones focales
help(focal)
promedio <- focal(r1, w = matrix(1/9, ncol = 3, nrow = 3))
promedio <- focal(r1, w = matrix(1, 3, 3), fun = "mean") # otra forma
print(as.matrix(r1))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 1 2 2
FALSE [2,] 1 1 2 2
FALSE [3,] 1 1 3 3print(as.matrix(promedio))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] NaN NaN NaN NaN
FALSE [2,] NaN 1.444444 1.888889 NaN
FALSE [3,] NaN NaN NaN NaNEl paquete raster permite también llevar a cabo operaciones zonales que consiste en calcular algún indice sobre un mapa, estratificando el cálculo con base en otro mapa. En el ejemplo a continuación, se suman el valor de las celdas del mapa r2 para cada categoría (valor) del mapa r1 (Figura 5). El resultado obtenido es una tabla.
# Operaciones zonales
print(as.matrix(r1))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 1 2 2
FALSE [2,] 1 1 2 2
FALSE [3,] 1 1 3 3print(as.matrix(r2))FALSE [,1] [,2] [,3] [,4]
FALSE [1,] 1 2 9 9
FALSE [2,] 2 NA 2 2
FALSE [3,] 9 9 9 1zonal_sum <- zonal(r2, r1, "sum") # map de valores, mapa de zonas
print(zonal_sum)FALSE zone sum
FALSE [1,] 1 23
FALSE [2,] 2 22
FALSE [3,] 3 106.2
En esta sección, vamos a aplicar algunos de los métodos anteriores para elaborar un mapa de áreas deforestadas basados en dos mapas de cubierta/uso del suelo de una región de Michoacán para 2004 y 2014 (Mas et al., 2017).
Para ello vamos a importar y rasterizar los mapas que se encuentran en formato shape utilizando el campo GRIDCODE en la tabla de atributos del mapa en formato vector y una resolución espacial de 100 m. La importación se realizó con la función st_read() de sf y por lo tanto se tuvo que convertir el objeto espacial sf a sp porque rasterize() solo acepta objetos de sp. En este caso, dependiendo de su computador, la rasterización puede tardar unos minutos. El proceso puede optimizarse utilizando la función gdal_rasterize del paquete gdalUtils.
# Importación de un shape con sf
library(sf)
cs2004 <- st_read("recursos-mx/cs2004.shp")FALSE Reading layer `cs2004' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\cs2004.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 4301 features and 5 fields
FALSE Geometry type: MULTIPOLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: 2460000 ymin: 780000 xmax: 2540000 ymax: 820000
FALSE Projected CRS: Lambert_Conformal_Conicplot(cs2004["Cl_abrev"], axes = T,cex.lab=0.8)
head(cs2004)FALSE Simple feature collection with 6 features and 5 fields
FALSE Geometry type: MULTIPOLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: 2463605 ymin: 780000 xmax: 2525585 ymax: 816174.4
FALSE Projected CRS: Lambert_Conformal_Conic
FALSE OBJECTID GRIDCODE clase area Cl_abrev
FALSE 1 396 1 Agricultura de riego 167.408719 Agr
FALSE 2 397 1 Agricultura de riego 27.766582 Agr
FALSE 3 398 1 Agricultura de riego 83.592112 Agr
FALSE 4 412 1 Agricultura de riego 2.912169 Agr
FALSE 5 418 1 Agricultura de riego 4.006450 Agr
FALSE 6 8596 2 Agricultura de temporal 1.966835 Agt
FALSE geometry
FALSE 1 MULTIPOLYGON (((2525547 780...
FALSE 2 MULTIPOLYGON (((2497969 780...
FALSE 3 MULTIPOLYGON (((2495969 780...
FALSE 4 MULTIPOLYGON (((2463803 791...
FALSE 5 MULTIPOLYGON (((2505625 793...
FALSE 6 MULTIPOLYGON (((2523434 816...st_bbox(cs2004)FALSE xmin ymin xmax ymax
FALSE 2460000 780000 2540000 820000st_crs(cs2004)FALSE Coordinate Reference System:
FALSE User input: Lambert_Conformal_Conic
FALSE wkt:
FALSE PROJCRS["Lambert_Conformal_Conic",
FALSE BASEGEOGCRS["WGS 84",
FALSE DATUM["World Geodetic System 1984",
FALSE ELLIPSOID["WGS 84",6378137,298.257223563,
FALSE LENGTHUNIT["metre",1]],
FALSE ID["EPSG",6326]],
FALSE PRIMEM["Greenwich",0,
FALSE ANGLEUNIT["Degree",0.0174532925199433]]],
FALSE CONVERSION["unnamed",
FALSE METHOD["Lambert Conic Conformal (2SP)",
FALSE ID["EPSG",9802]],
FALSE PARAMETER["Latitude of false origin",12,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8821]],
FALSE PARAMETER["Longitude of false origin",-102,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8822]],
FALSE PARAMETER["Latitude of 1st standard parallel",17.5,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8823]],
FALSE PARAMETER["Latitude of 2nd standard parallel",29.5,
FALSE ANGLEUNIT["Degree",0.0174532925199433],
FALSE ID["EPSG",8824]],
FALSE PARAMETER["Easting at false origin",2500000,
FALSE LENGTHUNIT["metre",1],
FALSE ID["EPSG",8826]],
FALSE PARAMETER["Northing at false origin",0,
FALSE LENGTHUNIT["metre",1],
FALSE ID["EPSG",8827]]],
FALSE CS[Cartesian,2],
FALSE AXIS["(E)",east,
FALSE ORDER[1],
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]],
FALSE AXIS["(N)",north,
FALSE ORDER[2],
FALSE LENGTHUNIT["metre",1,
FALSE ID["EPSG",9001]]]]## Rasteriza con paquete raster
resol <- 100
# Para rasterizar usando un raster "master"
extension <- extent(c(2460000, 2540000, 780000, 820000))
r <- raster(res = resol, ext = extension)
# Rasteriza con el campo GRIDCODE (as() para convertir a sp)
cs2004r <- rasterize(as(cs2004, Class = "Spatial"), r, "GRIDCODE")
plot(cs2004r, axes = F)
## Lo mismo con el mapa de 2014
cs2014 <- st_read("recursos-mx/cs2014.shp")FALSE Reading layer `cs2014' from data source
FALSE `G:\Mi unidad\2022-1\Estadistica Espacial\recursos-mx\cs2014.shp'
FALSE using driver `ESRI Shapefile'
FALSE Simple feature collection with 4347 features and 4 fields
FALSE Geometry type: MULTIPOLYGON
FALSE Dimension: XY
FALSE Bounding box: xmin: 2460000 ymin: 780000 xmax: 2540000 ymax: 820000
FALSE Projected CRS: Lambert_Conformal_Conic# plot(cs2014["GRIDCODE"], axes = T)
cs2014r <- rasterize(as(cs2014, Class = "Spatial"), r, "GRIDCODE")
plot(cs2014r, axes = F)
Elaboramos una tabla de correspondencia entre los valores de los campos GRIDCODE y clase realizando un resumen de la tabla de atributos con la función aggregate() que vimos en el capítulo 2 (sección 2.7). Note que este aggregate para el manejo de tablas pertenece al paquete stats y no es el mismo que el que realiza una agregación espacial de las celdas de un mapa raster. R sabe reconocer cual de los dos debe usar por el tipo de insumo en los argumentos de la función. Sin embargo, para evitar toda ambiguedad, se puede especificar el paquete utilizado como en aggregate(GRIDCODE~clase,FUN=min,data=cs2004,package=“stat”).
# Obtiene GRIS CODE / clase
# aggregate(GRIDCODE~clase,FUN=min,data=cs2004,package="stat")
aggregate(GRIDCODE~clase, FUN = min, data = cs2004) # equivalente a anteriorFALSE clase GRIDCODE
FALSE 1 Agricultura de riego 1
FALSE 2 Agricultura de temporal 2
FALSE 3 Asentamientos humanos 4
FALSE 4 Bosque de encino/veg primaria arbrea 6
FALSE 5 Bosque de encino/veg secundaria herbcea 7
FALSE 6 Bosque de pino/veg primaria 10
FALSE 7 Bosque de pino/veg secundaria 11
FALSE 8 Bosque mesofilo secundario 13
FALSE 9 Bosque Pino encino/veg primaria 14
FALSE 10 Bosque Pino encino/veg secundaria 15
FALSE 11 Cuerpos de agua 20
FALSE 12 Cultivo perenne 3
FALSE 13 Pastizal inducido pastizal cultivado 5
FALSE 14 Selva baja caducifolia/veg primaria 16
FALSE 15 Selva baja caducifolia/veg secundaria 17
FALSE 16 Selva mediana caducifolia/veg primaria 18
FALSE 17 Selva mediana caducifolia/veg secundaria 19
FALSE 18 Sin vegetacin aparente 23En seguida, vamos a reclasificar los mapas en dos categorías (forestal y no forestal) y combinar los mapas de 2004 y 2014 para generar un mapa de las áreas deforestadas durante este periodo. La reclasificación se base en una tabla (matriz con tres columnas) que indica los límites de los intervalos de reclasificación y el valor que tomarán las celdas.
# Reclasificación para elaborar un mapa binario forestal / no forestal
# 6-15 y 16-19 son clases forestales
# Matriz tabla para reclasificar
# por default el valor inferior del rango no está incluido
m <- c(0, 5, 0, 5, 15, 1, 15, 19, 1, 19, 24, 0)
rclmat <- matrix(m, ncol = 3, byrow = TRUE)
print(rclmat)FALSE [,1] [,2] [,3]
FALSE [1,] 0 5 0
FALSE [2,] 5 15 1
FALSE [3,] 15 19 1
FALSE [4,] 19 24 0F2004 <- reclassify(cs2004r, rclmat)
plot(F2004,axes=F)
F2014 <- reclassify(cs2014r, rclmat)
plot(F2014, axes = F)
## Algebra de mapas: Deforestación
def <- overlay(F2004, F2014, fun = function(x,y)
{ifelse( x == 1 & y == 0, 1, 0)})
plot(def, axes=F)
## Agrega pixels por paquetes de 10 x 10 pixels (suma)
def2 <- aggregate(def, fact=10, fun=sum)
plot(def2,axes=F)
# Salva el raster
writeRaster(def2, filename = "defor.tif", format = "GTiff", datatype = "INT1U",
overwrite = TRUE)Referencia
Mas, Jean. (2018). Análisis espacial con R: Usa R como un Sistema de Información Geográfica.