Capítulo 5 Graficar datos espaciales (SHP) en R
5.1 Objetivo de la practica
Comprender el procedimiento para importar archivos geoespaciales en formato shapefile en R e identificar los pasos básicos para representar mapas sencillos a partir de dichos datos.
5.2 Resultados de la práctica
Importar un archivo shp en R y generar un mapa básico que muestre la información contenida en dicho archivo, comprendiendo su utilidad como insumo para la generación ágil de mapas.
5.3 Instalación de librerias necesarias
Para este clase necesitaremos 1 librería que si aún no las has usado es importante que las instales previamente.
- librería
sf: permite manejar y analizar datos espaciales en R mediante el estándar simple features, facilitando la lectura, transformación y visualización de información geográfica.
5.4 Que datos usaremos en la clase
Para esta práctica se utilizarán dos datos de dos tipos diferentes. El primer insumo corresponde al límite del Distrito Metropolitano de Quito en formato shapefile (SHP), que permitirá la representación espacial del área de estudio. El segundo insumo es un archivo en formato CSV, que contiene las coordenadas (x, y) de 17 puntos: 9 estaciones de la Red Metropolitana de Monitoreo Atmosférico de Quito y 8 puntos de referencia obtenidos de la base ERA5. Estos datos se emplearán para ilustrar la integración de información vectorial y tabular en un entorno SIG desde R.
5.5 Empecemos
Antes de empezar cualquier proceso que requiera un sistema de referencia de coordenadas para trabajar, es necesario y muy importante definir este sistema que utilizaremos para nuestro análisis. Es crucial que todos los recursos y datos estén en el mismo sistema de referencia para asegurar la coherencia y precisión en el análisis. Si los datos provienen de diferentes fuentes y tienen sistemas de referencia distintos, es necesario transformarlos al sistema de referencia común. En esta ocación se empleará el sistema de coordenadas UTM.
En la línea de codigo que veremos a continuación, únicamente estamos almacenando en una variable el sistema de coordenadas con el que trabajaremos durante el proyecto; tenerlo como una variable nos ahorrará mucho tiempo de escritura cada vez que lo necesitemos.
utm17s <- "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs"Emplearemos la función shapefile(), de la librería sp, para cargar archivos tipo shapefile (SHP) y posteriormente con la función print() visualizaremos la información que contiene el shape, este paso es fundamental ya que aquí podremos comprobar el sistema de referencia de coordenadas en el que está el archivo.
DMQ<-st_read("./Documentos/DMQ_arc/limite_dmq.shp")
#> Reading layer `limite_dmq' from data source
#> `D:\PAO\2025B\Hidrologia_Aplicada\Hidrologia_aplicada\Documentos\DMQ_arc\limite_dmq.shp'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 450069.9 ymin: 9934660 xmax: 536741.5 ymax: 10028270
#> Projected CRS: quitow84
print(DMQ) # Nos mostrará la información del shapefile
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 450069.9 ymin: 9934660 xmax: 536741.5 ymax: 10028270
#> Projected CRS: quitow84
#> NOMBRE geometry
#> 1 limite del DMQ POLYGON ((504335.4 10027355...Observemos que el CRS (Sistema de referencia de coordenadas) no corresponde a UTM, por lo que será necesario cambiar de sistema al que hemos establecido como CRS del proyecto.
Al cargar un archivo vectorial con el paquete sf, este se convierte en un objeto de tipo sf (simple features). A estos objetos es necesario verificar y asignarles un sistema de referencia de coordenadas (CRS). Primero se debe declarar el CRS original (si no viene definido en el archivo) y, posteriormente, se puede transformar a otro sistema de coordenadas.
st_crs()permite consultar o asignar un sistema de coordenadas a un objeto sf.st_transform()se utiliza para cambiar el sistema de coordenadas de un objeto sf a otro CRS
En las líneas de código anteriores verificamos que el archivo tenga definido un SRC y además observamos que este no corresponde al sistema de referencia del proyecto. En esta ocación no será necesario utilizar la función st_crs(), pero es muy importante transformar el SRC.
Para imprimir el archivo con la geometría, usaremos dos funciones: plot() y st_geometry(). Recuerda que plot ya lo hemos usado varias veces y es el encargado de de graficas, por otra parte st_geometry extrae la geometría del objeto.
DMQ<- st_transform(DMQ, crs = utm17s)
print(DMQ)
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 728348 ymin: 9934644 xmax: 815024.9 ymax: 10028270
#> Projected CRS: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
#> NOMBRE geometry
#> 1 limite del DMQ POLYGON ((782609.7 10027358...
#Imprimir el archivo de datos espaciales
plot(st_geometry(DMQ), axes=TRUE, main = "Distrito Metropolitano de Quito")
Ahora carguemos el archivo txt que contiene la ubicación de 17 estaciones. Observa que usamos la función read.csv() para leer el archivo CSV donde los valores están separados por comas.
EST.XY <- read.csv("./Documentos/Estaciones.csv")
head(EST.XY)
#> ESTACION ALT LON LAT
#> 1 BEL 2797 779400 9980085
#> 2 CAR 2659 784166 9989120
#> 3 CEN 2820 777160 9975660
#> 4 COT 2752 778596 9988076
#> 5 CAM 2922 777171 9972341
#> 6 GUA 3047 772559 9963399st_a_sf() convierte un data frame regular (EST.XY) en un objeto de tipo sf al asignarle coordenadas espaciales.Además observa que asignamos al objeto un sistema de coordenadas proyectadas.
Una vez que todos nuestros datos estén en el mismo sistema de referencia, resulta sencillo generar gráficos que combinan tanto las estaciones como el contorno del Distrito Metropolitano de Quito (DMQ). La homogeneidad en el sistema de referencia facilita la superposición de las capas de datos.
Grafiquemos un mapa donde se visualice el polígono del Distrito Metropolitano de Quito y, sobre él, la ubicación y nombre de las estaciones.
Este código dibuja primero el límite del Distrito Metropolitano de Quito usando plot(st_geometry(DMQ)), añadiendo ejes y un título al mapa. Luego, con plot(st_geometry(EST.XY), add=TRUE) se colocan sobre el mapa los puntos que representan las estaciones, usando un símbolo circular sencillo. Finalmente, con text(st_geometry(EST.XY)) se añaden las etiquetas de cada estación justo encima de su punto, tomando los nombres de la columna ESTACION.
# Graficar el polígono del DMQ
plot(st_geometry(DMQ), axes = TRUE, main = "Distrito Metropolitano de Quito")
# Agregar las estaciones como puntos
plot(st_geometry(EST.XY), col = "black", pch = 21, add = TRUE)
# Agregar etiquetas a los puntos
text(st_coordinates(EST.XY), labels = EST.XY$ESTACION, pos = 3, cex = 0.7, col = "black")