Informe técnico
Simón Pedro Galeano Muñoz
Capítulo 3: Organización de los objetos espaciales en R
En principio hacen una introducción acerca de los objetos presentados en el capítulo, los cuales (claro está, objetos para representar infor mación espacial) son objetos vectoriales y objetos de tipo raster.
Dicho esto, el autor hace alusión a que se usará el paquete sf puesto que esté paulatinamente será el reemplazo del paquete sp, además, presenta razones por las cuales es recomendable usar el paquete sf.
Posteriormente, se trae a colación el hecho de que si bien en la práctica los datos espaciales son cargados y no construidos (por lo general), es necesario e importante entender como funcionan dichos datos y por tanto se procede con una exploración de algunos objetos de la librería sf.
Cobertura de puntos
Creación de puntos
#st_point se usa para un solo punto
P <- st_point(c(2, 5)) #punto en dos dimensiones (X, Y)
#se tiene un maximo de 4 dimensiones
print(class(P)) #viendo las clases a las que pertenece el objeto P## [1] "XY" "POINT" "sfg"str(P) #representacion como string del objeto P## 'XY' num [1:2] 2 5plot(P, axes=T, xlab="x", ylab="y",
pch=18, main="Representación del punto en el plano",
xlim = c(-10, 10), ylim=c(-10, 10)) 
#coordenadas en X e Y
Xs <- c(2, 4, 5); Ys <- c(5, 4, 8)
coords <- cbind(Xs, Ys) #pegando por columnas
print(coords) #mostrndo el resultado## Xs Ys
## [1,] 2 5
## [2,] 4 4
## [3,] 5 8MP <- st_multipoint(coords) #st_multipoint es para varios puntos
plot(MP, axes=T, xlab="x", ylab="y",
pch=18, main="Representación de los puntos en el plano",
xlim = c(-10, 10), ylim=c(-10, 10))
print(class(MP)) #analogo a lo que se hizo con P## [1] "XY" "MULTIPOINT" "sfg"print(MP) #representacion separada por pares ordenados del punto## MULTIPOINT ((2 5), (4 4), (5 8))#puntos en 3d
xyz <- cbind(coords, Zs=c(17, 22, 31))
print(xyz) #mostrando las ternas (puntos)## Xs Ys Zs
## [1,] 2 5 17
## [2,] 4 4 22
## [3,] 5 8 31MP3 <- st_multipoint(xyz) #convirtiendo las ternas en multipuntos
print(MP3)## MULTIPOINT Z ((2 5 17), (4 4 22), (5 8 31))Colleciones sfc y su manipulación
La idea de este apartado es mostrar como crear colecciones de puntos a partir de puntos individuales
#creando puntos
P1 <- st_point(c(2, 5)); P2 <- st_point(c(4, 4)); P3 <- st_point(c(5, 8))
geometria1 <- st_sfc(P1, P2, P3) #cranto la coleccion de puntos
#con el argumento crs se puede especificar proyeccion geograficaLuego, se busca construir un objeto más sofisticado como se muestra a continuación
#generando tabla de atributos
num <- 1:3
nombre <- c("Pozo", "Gasolinera", "Pozo")
#data frame con el id y la clase de lugar
tabpuntos <- data.frame(cbind(num, nombre))
#class(tabpuntos) da como salida data.frame
tabpuntos## num nombre
## 1 1 Pozo
## 2 2 Gasolinera
## 3 3 Pozo#creando un objeto de caracteristicas simples
SFP <- st_sf(tabpuntos, geometry=geometria1)
class(SFP) #viendo las clases a las que pertence el objeto## [1] "sf" "data.frame"print(SFP) #observando la combinacion de atributos y geometria## Simple feature collection with 3 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 2 ymin: 4 xmax: 5 ymax: 8
## CRS: NA
## num nombre geometry
## 1 1 Pozo POINT (2 5)
## 2 2 Gasolinera POINT (4 4)
## 3 3 Pozo POINT (5 8)plot(SFP, axes=T) #graficando los puntos por los atributos (id y nombre del sitio)
#st_geometry extrae la geometria del objeto sp y luego esta es graficada
plot(st_geometry(SFP), axes=T, xlab="x", ylab="y", main="Gráfico de la geometría", pch=20)
Nota: se puede filtrar usando la sintáxis usual de fancy indexing de R, además, es posible tener solo la parte de los atributos usando la función as.data.frame en el objeto sf.
as.data.frame(SFP) #extrayendo solo los atributos## num nombre geometry
## 1 1 Pozo POINT (2 5)
## 2 2 Gasolinera POINT (4 4)
## 3 3 Pozo POINT (5 8)Cobertura de líneas
La idea es unir pares ordenados \(\{(x_i, y_i)\}_{i = 1}^{n}\) a traves del segmento de recta que se genera enter dos puntos sucesivos de la colección previamente descrita, siempre y cuando estos segmentos de recta unidos no presenten ni bifurcaciones ni intersecciones, esto a traves de la función st_linestring.
X1s <- c(0, 3, 5, 8, 10)
Y1s <- c(0, 3, 4, 8, 10)
Coord1 <- cbind(X1s, Y1s)
L1 <- st_linestring(Coord1) #creando la primera union de segmentos (linea)
class(L1) #viendo las clases de dicho objeto## [1] "XY" "LINESTRING" "sfg"print(L1) #pares ordenados por los que pasa la linea## LINESTRING (0 0, 3 3, 5 4, 8 8, 10 10)X2s <- c(2, 1, 1)
Y2s <- c(2, 4, 5)
Coord2 <- cbind(X2s, Y2s)
L2 <- st_linestring(Coord2) #segunda linea
X3s <- c(8, 8)
Y3s <- c(8, 5)
Coord3 <- cbind(X3s, Y3s)
L3 <- st_linestring(Coord3)Estos objetos se pueden agrupar de una manera similar a como se agrupan los puntos, como se muestra a continuación.
L1L2L3 <- st_multilinestring(list(L1, L2, L3)) #note que se pasan en una lista y no separados por comas
plot(L1L2L3, axes=T) #lineas agrupadas
#se podria producir el mismo resultado haciendo plot(L1, axes=T); plot(L2, add = T); plot(L2, add = T)Asimismo, se pueden combinar geometrías con atributos no necesariamente geométricos
num <- 1:3
code <- c("t", "t", "p")
tipo <- c("Terraceria", "Terraceria", "Pavimentada")
tablineas <- data.frame(cbind(num, tipo, code))
tablineas #viendo los atributos## num tipo code
## 1 1 Terraceria t
## 2 2 Terraceria t
## 3 3 Pavimentada pgeometria2 <- st_sfc(L1, L2, L3) #uniendo las lineas en una coleccion
SFL <- st_sf(tablineas, geometry=geometria2) #objeto de caracteristicas simples
#graficando por atributos y coloreando por las clases dentro de los atributos
plot(SFL, axes=T)
De manera análoga como con los puntos
as.data.frame(SFL) #extrayendo la tabla de atributos## num tipo code geometry
## 1 1 Terraceria t LINESTRING (0 0, 3 3, 5 4, ...
## 2 2 Terraceria t LINESTRING (2 2, 1 4, 1 5)
## 3 3 Pavimentada p LINESTRING (8 8, 8 5)SFL[tipo=="Pavimentada", ] #filtrando con la sintaxis nativa de R## Simple feature collection with 1 feature and 3 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 8 ymin: 5 xmax: 8 ymax: 8
## CRS: NA
## num tipo code geometry
## 3 3 Pavimentada p LINESTRING (8 8, 8 5)plot(st_geometry(SFL[tipo=="Terraceria", ]), axes=T, col="red", xlab="x", ylab="y",
main = "Segmentos distinguidos por tipo")
plot(st_geometry(SFL[tipo=="Pavimentada", ]), add=T)
Cobertura de polígonos
Se crean polígonos usando sus vértices con la función st_polygon, dicho polígono debe ser cerrado. Si la secuencia se mueve en sentido horario, es para generar agujeros en otro polígono, mientras que si está en sentido contrario se indica que se refiere a la creación de un contorno externo.
X1 <- c(5, 10, 10, 6, 5); Y1 <- c(0, 0, 5, 5, 0)
c1 <- cbind(X1, Y1); print(c1) #mostrando un objeto de la clase matrix con los vertices## X1 Y1
## [1,] 5 0
## [2,] 10 0
## [3,] 10 5
## [4,] 6 5
## [5,] 5 0P1 <- st_polygon(list(c1)) #noe que se pasa como una lista
plot(P1, axes=T, xlab="x", ylab="y", main="Polígono")
#contorno del poligono
X2 <- c(0, 2, 6, 0, 0); Y2 <- c(4, 4, 10, 10, 4)
c2 <- cbind(X2, Y2)
#agujero al interior del poligono
X3 <- c(1, 1, 2, 2, 1); Y3 <- c(5, 6, 6, 5, 5)
c3 <- cbind(X3, Y3)
P2 <- st_polygon(list(c2, c3))
plot(P2, axes=T, xlab="x", ylab="y", main="Polígono hueco")
Un solo polígono por supuesto que no es hueco
plot(P4 <- st_polygon(list(c3[nrow(c3):1, ])),
axes=T, xlab="x", ylab="y", main="Polígono no hueco")
Otro ejemplo de polígono
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, xlab="x", ylab="y", main="Polígono de área urbana")
Ahora, se ve como se pueden asignar también atributos y clear un objeto de atributos simples.
ID <- 1:4; code <- c("F", "F", "U" , "A")
tipo <- c("Bosque", "Bosque", "Urbano", "Agricultura")
tabpol <- data.frame(cbind(ID, code, tipo))
#objeto sf
geometria3 <- st_sfc(P1, P2, P4, P5) #coleccion de atributos simples
SFPol <- st_sf(tabpol, geometry=geometria3)
plot(SFPol, axes=T) #mas de lo mismo 
Como se ha venido haciendo hasta ahora, se muestran algunos usos adicionales de el objeto creado con los polígonos y atributos.
class(SFPol) #clases del objeto## [1] "sf" "data.frame"print(SFPol) #caracteristicas del objeto## Simple feature collection with 4 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
## CRS: NA
## ID code tipo geometry
## 1 1 F Bosque POLYGON ((5 0, 10 0, 10 5, ...
## 2 2 F Bosque POLYGON ((0 4, 2 4, 6 10, 0...
## 3 3 U Urbano POLYGON ((1 5, 2 5, 2 6, 1 ...
## 4 4 A Agricultura POLYGON ((0 0, 5 0, 6 5, 10...summary(SFPol) #(nuevo) resumen del objeto## ID code tipo geometry
## Length:4 Length:4 Length:4 POLYGON:4
## Class :character Class :character Class :character epsg:NA:0
## Mode :character Mode :character Mode :characteras.data.frame(SFPol) #extrayendo los atributos## ID code tipo geometry
## 1 1 F Bosque POLYGON ((5 0, 10 0, 10 5, ...
## 2 2 F Bosque POLYGON ((0 4, 2 4, 6 10, 0...
## 3 3 U Urbano POLYGON ((1 5, 2 5, 2 6, 1 ...
## 4 4 A Agricultura POLYGON ((0 0, 5 0, 6 5, 10...print(Bosque <- SFPol[tipo=="Bosque", ])## Simple feature collection with 2 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
## CRS: NA
## ID code tipo geometry
## 1 1 F Bosque POLYGON ((5 0, 10 0, 10 5, ...
## 2 2 F Bosque POLYGON ((0 4, 2 4, 6 10, 0...#graficando la geometria del objeto
plot(st_geometry(SFPol), axes=T, xlab="x", ylab="y", main="Bosque")
plot(st_geometry(Bosque), col="green", add=T)
Diferentes tipos de geometría
Así como se ha venido creando colecciones de objetos del mismo tipo, también se pueden hacer colecciones de objetos con distintos tipos de geometría.
Detodo <- st_sfc(P1, P2, L1, L2)
plot(Detodo, border="red", axes=T, col="black", lwd=3,
xlab="x", ylab="y", main="Combinando lineas y polígonos")
Datos raster
También se pueden manejar datos raster a través del paquete raster, obteniendo el objeto a partir de una matriz con los valores de cada celda de dicho arreglo.
library(raster)## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectm <- matrix(c(1, 2, 3, 4, 2, NA, 2, 2, 3, 3, 3, 1), ncol=4, nrow=3, byrow=T)
m## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 2 NA 2 2
## [3,] 3 3 3 1Se muestra como es la creación del objeto raster
r <- raster(m)
extent(r) <- extent(c(0, 4, 0, 3)) #dimensiones del raster
class(r) #clase del objeto raster## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"print(r) #representacion del objeto raster## class : RasterLayer
## dimensions : 3, 4, 12 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0, 4, 0, 3 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : layer
## values : 1, 4 (min, max)plot(r, xlab="x", ylab="y", main="Objeto Raster")
Capítulo 4: Importación/exportación de datos espaciales
Importación de archivos shape
Los archivos shape son aquellos que tiene por extensión shapefile, el cual se encuentra compuesto por lo menos por tres archivos con las siguientes extensiones y contenidos:
- .shp: son las entidades geométricas de los objetos.
- .shx: se encarga de almacenar la indexación de las entidades geométricas.
- .dbf: almacena los atributos de los objetos.
Luego de esto se hace mención de la existencia de múltiples maneras de realizar la importación de estos archivos en R, sin embargo como viene siendo de costumbre a lo largo de la lectura, se ilustra como realizar dicha tarea con el paquete sf:
col <- st_read("COL_adm\\COL_adm1.shp") #leyendo .shp## Reading layer `COL_adm1' from data source
## `C:\Users\sigal\OneDrive\Desktop\SemestreVIII\EstadisticaEspacial\AnalisisEspacialR\COL_adm\COL_adm1.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.84153 ymin: -4.228429 xmax: -66.87033 ymax: 15.91247
## Geodetic CRS: WGS 84#verificando que es un sf
class(col)## [1] "sf" "data.frame"summary(col) #resumen de la base de datos## ID_0 ISO NAME_0 ID_1
## Min. :53 Length:32 Length:32 Min. : 1.00
## 1st Qu.:53 Class :character Class :character 1st Qu.: 8.75
## Median :53 Mode :character Mode :character Median :16.50
## Mean :53 Mean :16.50
## 3rd Qu.:53 3rd Qu.:24.25
## Max. :53 Max. :32.00
## NAME_1 TYPE_1 ENGTYPE_1 NL_NAME_1
## Length:32 Length:32 Length:32 Length:32
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## VARNAME_1 geometry
## Length:32 MULTIPOLYGON :32
## Class :character epsg:4326 : 0
## Mode :character +proj=long...: 0
##
##
## plot(st_geometry(col), axes=T) #solo se grafica la geometria porque el .shp tiene bsatantes columnas
Ahora, se observan al gunas propiedades del objeto shape
st_geometry(col) #geometria del objeto## Geometry set for 32 features
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.84153 ymin: -4.228429 xmax: -66.87033 ymax: 15.91247
## Geodetic CRS: WGS 84
## First 5 geometries:## MULTIPOLYGON (((-69.43138 -1.078474, -69.42712 ...## MULTIPOLYGON (((-76.92486 8.212361, -76.92486 8...## MULTIPOLYGON (((-70.6878 7.09091, -70.66744 7.0...## MULTIPOLYGON (((-74.84764 11.10792, -74.84764 1...## MULTIPOLYGON (((-75.67486 10.16125, -75.67486 1...st_crs(col) #sistema de coordenadas de referencia## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## 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",4326]]También se puede añadir el sistema de coordenadas en caso tal de que este no se encuentre, esto a través de una asignación, adicionalmente, se puede verificar la validez de los mapas usando la función st_is_valid.
Importando archivos de otros formatos distintos
Se presentan algunas clases diferenctes de archvos que pueden ser leídos con la función st_read diferentes al clásico .shp, con información adicional como si se puede escribir, copiar, si son vectores, etc.
head(st_drivers()[, -2], n=10)## name write copy is_raster is_vector vsi
## ESRIC ESRIC FALSE FALSE TRUE TRUE TRUE
## FITS FITS TRUE FALSE TRUE TRUE FALSE
## PCIDSK PCIDSK TRUE FALSE TRUE TRUE TRUE
## netCDF netCDF TRUE TRUE TRUE TRUE FALSE
## PDS4 PDS4 TRUE TRUE TRUE TRUE TRUE
## VICAR VICAR TRUE TRUE TRUE TRUE TRUE
## JP2OpenJPEG JP2OpenJPEG FALSE TRUE TRUE TRUE TRUE
## JPEG2000 JPEG2000 FALSE TRUE TRUE TRUE TRUE
## PDF PDF TRUE TRUE TRUE TRUE FALSE
## MBTiles MBTiles TRUE TRUE TRUE TRUE TRUEExportando a otros formatos
Se puede usar la función st_write y especificar el forato de escritura para un objeto ya existente que sea manipulable con la librería sf. Además, también se pueden exportar datos raster.
#guardando un objeto en formato shape
para_guardar <- st_sfc(P1, P2)
st_write(para_guardar, "detodo.shp", delete_layer=T)Adicionalmente se muestra el guardado de un archivo en un formato distinto al shape
#guardando en formato gpkg
st_write(para_guardar, "noraster.shp", delete_layer=T)Y para dar cierre, la librería sf permite convertir objetos sp sp a objectos sf
library(sp)
Detodo_sp <- as(Detodo, Class = "spatial") #objeto sp
Detodo_sf <- st_as_sf(Detodo_sp, "sf") #de objeto sp a sfImportación y exportación de datos raster
Como se ha venido haciendo mención, los objetos raster son de gran utilidad en el análisis de datos espaciales y en este apartado se muestra como importar y exportar dichos objetos.
dem <- raster("Map250k-full-300dpi.tif") #leyendo el objeto
class(dem) #clase del objeto raster## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"Se procede a visualizar el objeto
plot(dem)
Luego, para ver qué formatos están disponibles se puede usar la función writeFormats y con alguno de esos formatos guardar el raster.
writeRaster(dem, filename = "raster_propio.lan", format="LAN",
datatype="INT2S") #guardando el raster en LANCapítulo 5: Operaciones básicas de SIG (vector)
En principio el autor realiza la lectura de unos archivos, sin embargo eso se omite ya que dichos ficheros son previamente creados en este código y resulta redundante cargarlos de nuevo ya que no representan un espacio significativo en memoria.
Se pueden verificar propiedades deseables de los objetos sf, como por ejemplo qué geometrías son vaálidas, si las líneas son simples o no, el área de los polígonos, la medida de las líneas, entre otros.
#validez de poligonos
st_is_valid(SFPol)## [1] TRUE TRUE TRUE TRUE#lineas simples
st_is_simple(SFL)## [1] TRUE TRUE TRUE#area de poligonos
st_area(SFPol)## [1] 22.5 23.0 1.0 53.5#longitud de lineas
st_length(SFL)## [1] 14.307136 3.236068 3.000000#graficando los datos (NO ACORDE AL LIBRO!!!!)
plot(st_geometry(SFL), lty=3, lwd=2, col="red", axes=T)
plot(st_geometry(SFP), add=T)
plot(st_geometry(SFP) + c(0.5, 0), add=T, pch=c(49, 50, 51), cex=1)
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)
plot(SFPol["ID"], col = c("Dark Green", "green", "grey", "white"),
add=T)
Se puede evaluar si objetos se intersectan de este modo
#verificando intersecciones entre las lineas y poligonos
st_intersects(SFL, SFPol, sparse = F)## [,1] [,2] [,3] [,4]
## [1,] FALSE FALSE FALSE TRUE
## [2,] FALSE TRUE TRUE TRUE
## [3,] TRUE FALSE FALSE TRUEEsto retorna una matriz con 3 filas (cantidad de líneas) y 4 columnas (cantidad de polígonos), donde la entrada \(i, j\) de la matriz es TRUE si la línea \(i\) y el el polígono \(j\) se intersectan, además, tiene la opción de presentarse la matriz como una matriz sparse, ya que al ser booleana, resulta útil solo presentar las entradas que son verdaderas por fila ya que la otra información no es aportante a la hora de operar con dichas matrices.
También, se puede encontrar la distancia más corta entre elementos de dos objetos
st_distance(SFP, SFPol) #distancia entre poligonos y lineas## [,1] [,2] [,3] [,4]
## [1,] 3.922323 0.0000000 0.000000 0.5547002
## [2,] 1.765045 1.6641006 2.236068 0.0000000
## [3,] 3.162278 0.2773501 3.605551 0.0000000Nuevamente, la entrada \(i, j\) de la matriz denota la distancia más corta entre la línea \(i\) y el el polígono \(j\).
Finalmente, se pueden intersectar objetos como se muestra a continuación.
#interseccion de poligonos y puntos
SFPxSFPol <- st_intersection(SFP, SFPol)
print(SFPxSFPol)## Simple feature collection with 4 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 2 ymin: 4 xmax: 5 ymax: 8
## CRS: NA
## num nombre ID code tipo geometry
## 1 1 Pozo 2 F Bosque POINT (2 5)
## 1.1 1 Pozo 3 U Urbano POINT (2 5)
## 2 2 Gasolinera 4 A Agricultura POINT (4 4)
## 3 3 Pozo 4 A Agricultura POINT (5 8)#interseccionentre lineas y poligonos
SFLxSFPol <- st_intersection(SFL, SFPol)
print(SFLxSFPol)## Simple feature collection with 6 features and 6 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 0 ymin: 0 xmax: 10 ymax: 10
## CRS: NA
## num tipo code ID code.1 tipo.1 geometry
## 3 3 Pavimentada p 1 F Bosque POINT (8 5)
## 2 2 Terraceria t 2 F Bosque LINESTRING (1 4, 1 5)
## 2.1 2 Terraceria t 3 U Urbano POINT (1 5)
## 1 1 Terraceria t 4 A Agricultura LINESTRING (0 0, 3 3, 5 4, ...
## 2.2 2 Terraceria t 4 A Agricultura LINESTRING (2 2, 1 4)
## 3.1 3 Pavimentada p 4 A Agricultura LINESTRING (8 8, 8 5)Observando los resultados de las intersecciones
plot(SFLxSFPol["tipo.1"], lwd=3)
Se pueden extraer elementos de la intersección y hacer operaciones sobre estos.
SFLxSFPol_l <- SFLxSFPol %>%
st_collection_extract(type="LINESTRING") %>%
st_length()Análisis espacial en formato vector
mx <- st_read(".//recursos-mx/Entidades_latlong.shp")## Reading layer `Entidades_latlong' from data source
## `C:\Users\sigal\OneDrive\Desktop\SemestreVIII\EstadisticaEspacial\AnalisisEspacialR\recursos-mx\Entidades_latlong.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -118.4076 ymin: 14.5321 xmax: -86.71041 ymax: 32.71865
## Geodetic CRS: WGS 84plot(st_geometry(mx), axes=T)
Se pueden usar funciones como sts_crs para ver el sistema de coordenadas de referencia, también la clásica función summary para mostrar un resuen del objeto sf.
carr <- st_read(".//recursos-mx//carretera.shp")## Reading layer `carretera' from data source
## `C:\Users\sigal\OneDrive\Desktop\SemestreVIII\EstadisticaEspacial\AnalisisEspacialR\recursos-mx\carretera.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 12424 features and 4 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 1071375 ymin: 324561 xmax: 4076832 ymax: 2349296
## Projected CRS: Lambert_Conformal_Conicst_crs(carr) #sistema de coordenadas de referencia## Coordinate Reference System:
## User input: Lambert_Conformal_Conic
## wkt:
## PROJCRS["Lambert_Conformal_Conic",
## BASEGEOGCRS["GCS_GRS 1980(IUGG, 1980)",
## DATUM["D_unknown",
## ELLIPSOID["GRS80",6378137,298.257222101,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",12,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-102,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",17.5,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",29.5,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",2500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]summary(carr)## ADMINISTRA ENTIDAD NO_CARRILE SHAPE_len
## Length:12424 Length:12424 Length:12424 Min. : 2.23
## Class :character Class :character Class :character 1st Qu.: 1346.58
## Mode :character Mode :character Mode :character Median : 3780.05
## Mean : 6843.31
## 3rd Qu.: 8942.17
## Max. :118873.89
## geometry
## LINESTRING :12424
## epsg:NA : 0
## +proj=lcc ...: 0
##
##
## El autor propone hacer un análisis en conjunto de el territorio nacional con su sistema víal, en particular se usan las carreteras como objetivo principal del análisis.
mx_lcc <- st_transform(mx, st_crs(carr)) #transformando el sistema de coordenadas de referencia
st_crs(mx_lcc) #el mismo sistema que las carreteras## Coordinate Reference System:
## User input: Lambert_Conformal_Conic
## wkt:
## PROJCRS["Lambert_Conformal_Conic",
## BASEGEOGCRS["GCS_GRS 1980(IUGG, 1980)",
## DATUM["D_unknown",
## ELLIPSOID["GRS80",6378137,298.257222101,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",12,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-102,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",17.5,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",29.5,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",2500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]mx_lcc$AreaEdo <- st_area(mx_lcc) / 1000000 #se divide entre 1000^2 para obtener km^2head(cbind(Estado=mx_lcc$NOM_ENT, Area=round(mx_lcc$AreaEdo, 4))) %>%
kable(col.names = c("Estado", "Áream en km^2"))| Estado | Áream en km^2 |
|---|---|
| Distrito Federal | 1486.481 |
| Guerrero | 63565.1129 |
| Mexico | 22226.6434 |
| Morelos | 4859.432 |
| Sinaloa | 56802.9736 |
| Baja California | 73535.9098 |
Se procede a ver en conjunto los dos espacios geográficos que fueron los portagonistas en el análisis que desarrola el autor.
plot(st_geometry(carr), col="blue", axes=T)
plot(st_geometry(mx_lcc), border="red", add=T)
Se realiza un filtro de las carreteras federales en cada estado del sureste mexicano.
carrFed <- carr[carr$ADMINISTRA == "Federal", ] #seleccionando las carreteras federales
#obteniendo el sureste mexicano
sureste <- mx_lcc[mx_lcc$NOM_ENT %in% c("Tabasco", "Campeche", "Quintana Roo", "Yucatan"), ]
plot(st_geometry(sureste))
Luego se intersecta el mapa de las carreteras federales con el de los estados del sureste
carrFedXedos <- st_intersection(carrFed, sureste)
plot(st_geometry(carrFedXedos), col="blue", axes=T)
plot(st_geometry(sureste), add=T, border="red")
Ahora se procede a calcular la longitud de las carreteras.
long <- carrFedXedos %>%
st_collection_extract(type="LINESTRING") %>% #convirtiendo las multilines en lines
st_length() / 1000 #distancia en kilometros
carrFedXedos$long <- long #asignando la longitud a la interseccion
head(long)## Units: [m]
## [1] 1.219532 1.133280 19.162872 14.348592 1.500849 4.092854Finalmente se hace un resumen de la longitud por estado.
aggregate(long ~ NOM_ENT, FUN=sum, data = carrFedXedos) %>%
kable(col.names = c("Estado", "Longitud [m]"))| Estado | Longitud [m] |
|---|---|
| Campeche | 1056.2609 [m] |
| Quintana Roo | 876.6589 [m] |
| Tabasco | 583.5174 [m] |
| Yucatan | 1117.8214 [m] |
Luego de esto, se pretende calcular la proporción del territorio que se encuentra a menos de 20km de una carretera federal.
buf <- st_buffer(carrFed, dist=20000) #buffer de 20km
plot(st_geometry(buf))
Se realiza la intersección correspondiente y el cómputo de la proporción.
carrbufXedos <- st_intersection(buf, sureste) #intersectando el buffer y los estados del sureste
carrbufXedos$area <- st_area(carrbufXedos) / 1000000 #area de la interseccion en km^2
Suma_area <- aggregate(area ~ NOM_ENT, FUN=sum, data=carrbufXedos)
sureste <- merge(sureste, Suma_area) #combinando los objetos por clave
sureste$prop <- sureste$area/sureste$AreaEdoPara dar fin a este capítulo, se tiene la creación de un índice como se muestra a continuación.
tab_pop <- read.csv(".//recursos-mx//Pop2010_INEGI.csv") #leyendo los datos de la poblacion
head(tab_pop) %>%
kable()| NumEdo | Estado | Poptotal | Hombre | Mujer |
|---|---|---|---|---|
| 1 | Aguascalientes | 1184996 | 576638 | 608358 |
| 2 | Baja California | 3155070 | 1591610 | 1563460 |
| 3 | Baja California Sur | 637026 | 325433 | 311593 |
| 4 | Campeche | 822441 | 407721 | 41472 |
| 5 | Coahuila de Zaragoza | 2748391 | 1364197 | 1384194 |
| 6 | Colima | 650555 | 32279 | 327765 |
#combinando los datos de la poblacion con los geoespaciales
tab_pop <- tab_pop %>%
rename(CveEdo = NumEdo)
mx_lcc <- merge(mx_lcc, tab_pop)
#calculando la densidad poblacional
mx_lcc$dens <- mx_lcc$Poptotal/mx_lcc$AreaEdo
#calculando el indice de masculinidad
mx_lcc$IM <- 100*mx_lcc$Hombre/mx_lcc$Mujeraggregate(IM ~ NOM_ENT, FUN=function(x) x, data=mx_lcc) %>%
arrange(desc(IM)) %>%
head() %>%
kable(col.names = c("Estado", "Índice de masculinidad")) | Estado | Índice de masculinidad |
|---|---|
| Campeche | 983.12355 |
| Baja California Sur | 104.44169 |
| Baja California | 101.80049 |
| Sonora | 101.26573 |
| Nayarit | 99.45494 |
| Nuevo Leon | 99.43907 |
aggregate(dens ~ NOM_ENT, FUN=function(x) x, data=mx_lcc) %>%
arrange(desc(dens)) %>%
head() %>%
kable(col.names = c("Estado", "Densidad poblacional"))| Estado | Densidad poblacional |
|---|---|
| Distrito Federal | 5954.3850 [1/m^2] |
| Mexico | 682.7779 [1/m^2] |
| Morelos | 365.7273 [1/m^2] |
| Tlaxcala | 294.3981 [1/m^2] |
| Aguascalientes | 213.1794 [1/m^2] |
| Guanajuato | 180.8308 [1/m^2] |
Para dar culmen, se guarda el archivo como .shp como se muestra a continuación.
st_write(mx_lcc, "mx_lcc.shp", delete_layer=T)Capítulo 6: Operaciones básicas de SIG (raster)
En este capítulo el objetivo del autor es ilustrar las diferentes operaciones que se pueden realizar con datos raster, iniciando con ejemplos sencillos para luego trabajar con datos elaborados.
Algunas operaciones de análisis espacial
En aras de ilustrar las operaciones básicas, se utilizan matrices no muy complejas para dicho propósito.
#Creacion de matrices para trabajar
m1 <- matrix(c(1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 3, 3),
ncol = 4, nrow = 3, byrow = T)
m2 <- matrix(c(1, 2, 3, 4, 2, NA, 2, 2, 3, 3, 3, 1),
ncol = 4, nrow = 3, byrow = T)
print(m1)## [,1] [,2] [,3] [,4]
## [1,] 1 1 2 2
## [2,] 1 1 2 2
## [3,] 1 1 3 3print(m2)## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 2 NA 2 2
## [3,] 3 3 3 1r1 <- raster(m1); r2 <- raster(m2)
extent(r1) <- extent(r2) <- extent(c(0,4,0,3))
#Creando tablas de colores y almacenandolas en una pila
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)
Se puede obtener cierta información acerca del objeto raster.
#Informacion del raster
res(r1) #resolucion del raster## [1] 1 1dim(r12) #tamano de la resolucion del raster## [1] 3 4 2xmax(r1) # valores de interes del raster, existen también xmax, ymin y ymax## [1] 4cellStats(r1, "sum") # sumando los pixeles## [1] 20cellStats(r1, "sd") # obteniendo la desviacion de los pixeles## [1] 0.7784989Las operaciones básicas se pueden hacer de manera natural, en caso de querer realizar operaciones definidas por el usuario, se puede usar la función overlay, que recibe una función con tantos parámetros como objetos raster se introduzcan como argumento a overlay.
#Se hacen operaciones elemento a elemento como en arreglos
sum1 <- r1 + 2; print(as.matrix(sum1))## [,1] [,2] [,3] [,4]
## [1,] 3 3 4 4
## [2,] 3 3 4 4
## [3,] 3 3 5 5sum2 <- r1 + 2*r2; print(as.matrix(sum2))## [,1] [,2] [,3] [,4]
## [1,] 3 5 8 10
## [2,] 5 NA 6 6
## [3,] 7 7 9 5sum3 <- overlay(r1, r2, sum1, fun=function(x, y, z){ x + 2*y -z} )
print(as.matrix(sum3))## [,1] [,2] [,3] [,4]
## [1,] 0 2 4 6
## [2,] 2 NA 2 2
## [3,] 4 4 4 0Ahora, se observan operaciones más orientadas hacia los SIG.
#Operaciones de sistema de informacion geografica
#Mosaico (merge) y recorte (crop)
# Extent de r1 (xmin, xmin, ymin, ymax)
#notar que los
extent(r1)## class : Extent
## xmin : 0
## xmax : 4
## ymin : 0
## ymax : 3m3 <- matrix(rep(5, 12),ncol=4,nrow=3,byrow=TRUE)
print(m3)## [,1] [,2] [,3] [,4]
## [1,] 5 5 5 5
## [2,] 5 5 5 5
## [3,] 5 5 5 5r3 <- raster(m3) #convirtiendo m3 en objeto tipo raster
# extent(xmin, xmax, ymin, ymax)
extent(r3) <- extent(c(0,4,3,6)) # Está al norte de r1
mosaico <- merge(r1,r3) #combinando objetos raster (pegando arriba)
extent(mosaico)## class : Extent
## xmin : 0
## xmax : 4
## ymin : 0
## ymax : 6print(as.matrix(mosaico))## [,1] [,2] [,3] [,4]
## [1,] 5 5 5 5
## [2,] 5 5 5 5
## [3,] 5 5 5 5
## [4,] 1 1 2 2
## [5,] 1 1 2 2
## [6,] 1 1 3 3extent_corte <- extent(c(1,3,2,4))
corte <- crop(mosaico,extent_corte) #recortando raster
print(as.matrix(corte))## [,1] [,2]
## [1,] 5 5
## [2,] 1 2Se puede hacer reclasificación en el objeto raster.
# Reclasificaciones
print(as.matrix(r2))## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 2 NA 2 2
## [3,] 3 3 3 1r2[r2 > 2] <- 6
print(as.matrix(r2))## [,1] [,2] [,3] [,4]
## [1,] 1 2 6 6
## [2,] 2 NA 2 2
## [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)## [,1] [,2] [,3]
## [1,] 0 2 0
## [2,] 2 5 1reclas <- reclassify(mosaico, tabla_reclas)
print(as.matrix(reclas))## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 1
## [2,] 1 1 1 1
## [3,] 1 1 1 1
## [4,] 0 0 0 0
## [5,] 0 0 0 0
## [6,] 0 0 1 1# Substitución de valores
print(as.matrix(r2))## [,1] [,2] [,3] [,4]
## [1,] 1 2 6 6
## [2,] 2 NA 2 2
## [3,] 6 6 6 1r2[r2 == 6] <- 9
print(as.matrix(r2))## [,1] [,2] [,3] [,4]
## [1,] 1 2 9 9
## [2,] 2 NA 2 2
## [3,] 9 9 9 1# con tabla
tab_subs <- data.frame(id=c(1,2,3),v=c(1,56,84))
print(tab_subs)## id v
## 1 1 1
## 2 2 56
## 3 3 84sub <- subs(r1, tab_subs)
print(as.matrix(r1))## [,1] [,2] [,3] [,4]
## [1,] 1 1 2 2
## [2,] 1 1 2 2
## [3,] 1 1 3 3print(as.matrix(sub))## [,1] [,2] [,3] [,4]
## [1,] 1 1 56 56
## [2,] 1 1 56 56
## [3,] 1 1 84 84Incluso se puede por ejemplo agrupar celdas y aplicar funciones para generar objetos de dimensión menor.
# Remuestreos
# Agregación espacial
agreg <- aggregate(mosaico, fact=2, fun="sum")
print(as.matrix(mosaico));print(as.matrix(agreg))## [,1] [,2] [,3] [,4]
## [1,] 5 5 5 5
## [2,] 5 5 5 5
## [3,] 5 5 5 5
## [4,] 1 1 2 2
## [5,] 1 1 2 2
## [6,] 1 1 3 3## [,1] [,2]
## [1,] 20 20
## [2,] 12 14
## [3,] 4 10agreg2 <- aggregate(mosaico, fact=2, fun=modal,na.rm = TRUE)
print(as.matrix(agreg2))## [,1] [,2]
## [1,] 5 5
## [2,] 1 5
## [3,] 1 2# 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)); print(as.matrix(resb))## [,1] [,2] [,3] [,4]
## [1,] 5 5 5 5
## [2,] 5 5 5 5
## [3,] 1 1 2 2
## [4,] 1 1 2 2
## [5,] 1 1 3 3## [,1] [,2] [,3] [,4]
## [1,] 5 5 5.0 5.0
## [2,] 5 5 5.0 5.0
## [3,] 3 3 3.5 3.5
## [4,] 1 1 2.0 2.0
## [5,] 1 1 2.9 2.9Realizando convoluciones.
# 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))## [,1] [,2] [,3] [,4]
## [1,] 1 1 2 2
## [2,] 1 1 2 2
## [3,] 1 1 3 3print(as.matrix(promedio))## [,1] [,2] [,3] [,4]
## [1,] NaN NaN NaN NaN
## [2,] NaN 1.444444 1.888889 NaN
## [3,] NaN NaN NaN NaNSe permite definir zonas y allí operar sobre dichas zonas.
# Operaciones zonales
print(as.matrix(r1)); print(as.matrix(r2))## [,1] [,2] [,3] [,4]
## [1,] 1 1 2 2
## [2,] 1 1 2 2
## [3,] 1 1 3 3## [,1] [,2] [,3] [,4]
## [1,] 1 2 9 9
## [2,] 2 NA 2 2
## [3,] 9 9 9 1zonal_sum <- zonal(r2,r1,"sum") # map de valores, mapa de zonas
print(zonal_sum)## zone sum
## [1,] 1 23
## [2,] 2 22
## [3,] 3 10Análisis espacial en formato raster
Se procede a la conversión de un archivo shape a un objeto raster, para posteriormente hacer análisis con dicho objeto.
#Importando un archivo shape
cs2004 <- st_read(".//recursos-mx//cs2004.shp")## Reading layer `cs2004' from data source
## `C:\Users\sigal\OneDrive\Desktop\SemestreVIII\EstadisticaEspacial\AnalisisEspacialR\recursos-mx\cs2004.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 4301 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2460000 ymin: 780000 xmax: 2540000 ymax: 820000
## Projected CRS: Lambert_Conformal_Conicplot(cs2004["Cl_abrev"], axes = T,cex.lab=0.8)
head(cs2004) #encabezado ## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2463605 ymin: 780000 xmax: 2525585 ymax: 816174.4
## Projected CRS: Lambert_Conformal_Conic
## OBJECTID GRIDCODE clase area Cl_abrev
## 1 396 1 Agricultura de riego 167.408719 Agr
## 2 397 1 Agricultura de riego 27.766582 Agr
## 3 398 1 Agricultura de riego 83.592112 Agr
## 4 412 1 Agricultura de riego 2.912169 Agr
## 5 418 1 Agricultura de riego 4.006450 Agr
## 6 8596 2 Agricultura de temporal 1.966835 Agt
## geometry
## 1 MULTIPOLYGON (((2525547 780...
## 2 MULTIPOLYGON (((2497969 780...
## 3 MULTIPOLYGON (((2495969 780...
## 4 MULTIPOLYGON (((2463803 791...
## 5 MULTIPOLYGON (((2505625 793...
## 6 MULTIPOLYGON (((2523434 816...st_bbox(cs2004) #rango de x e y del objeto shape## xmin ymin xmax ymax
## 2460000 780000 2540000 820000st_crs(cs2004) #sistema de coordenadas de referencia## Coordinate Reference System:
## User input: Lambert_Conformal_Conic
## wkt:
## PROJCRS["Lambert_Conformal_Conic",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",12,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-102,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",17.5,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",29.5,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",2500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]#Rasterizando con el paquete raster
resol <- 100 #resolucion del raster
# Para rasterizar usando un raster "master"
extension <- extent(c(2460000, 2540000,780000, 820000))
r <- raster(res=resol, ext=extension) #objeto raster vacio
# 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")## Reading layer `cs2014' from data source
## `C:\Users\sigal\OneDrive\Desktop\SemestreVIII\EstadisticaEspacial\AnalisisEspacialR\recursos-mx\cs2014.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 4347 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2460000 ymin: 780000 xmax: 2540000 ymax: 820000
## Projected CRS: Lambert_Conformal_Conic# plot(cs2014["GRIDCODE"], axes = T)
cs2014r <- rasterize(as(cs2014,Class = "Spatial"), r, "GRIDCODE")
plot(cs2014r,axes=F)
Resumiendo información del raster.
# Obtiene GRIDCODE / clase
aggregate(GRIDCODE~clase,FUN=min,data=cs2004) # equivalente a anterior## clase GRIDCODE
## 1 Agricultura de riego 1
## 2 Agricultura de temporal 2
## 3 Asentamientos humanos 4
## 4 Bosque de encino/veg primaria arbrea 6
## 5 Bosque de encino/veg secundaria herbcea 7
## 6 Bosque de pino/veg primaria 10
## 7 Bosque de pino/veg secundaria 11
## 8 Bosque mesofilo secundario 13
## 9 Bosque Pino encino/veg primaria 14
## 10 Bosque Pino encino/veg secundaria 15
## 11 Cuerpos de agua 20
## 12 Cultivo perenne 3
## 13 Pastizal inducido pastizal cultivado 5
## 14 Selva baja caducifolia/veg primaria 16
## 15 Selva baja caducifolia/veg secundaria 17
## 16 Selva mediana caducifolia/veg primaria 18
## 17 Selva mediana caducifolia/veg secundaria 19
## 18 Sin vegetacin aparente 23Finalmente se elabora un mapa binario forestal y no forestal.
# 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)## [,1] [,2] [,3]
## [1,] 0 5 0
## [2,] 5 15 1
## [3,] 15 19 1
## [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)