Proceso de interpolacion espacial de variables hidroclimatológicas mediante IDW
Cristian David Vargas
3/04/2022
R Markdown
En este Rpub se desarrolla un ejemplo para una interpolación espacial de precipitación mediante el método IDW en la cuenca media del Río Sinú en Colombia.
Para poder desarrollar la interpolación se requiere de los siguientes insumos:
-Shape de la cuenca.
-Shape de las estaciones que miden precipitación cercanas a la cuenca media del Río Sinú.
-Una malla de puntos equidistantes en formato shape.
-Información de los valores de precipitación en una escala de tiempo determinada.
Puedes descargar los insumos para desarrollar este ejercicio en el siguiente link:
https://drive.google.com/drive/folders/1gNRb1bWyoVih7B4g_FbuxZ2nX6JdMdyp?usp=sharing
Es importante que todos los shapes que se utilicen para la interpolación se encuentren en un mismo sistema de coordenadas planas.
En este ejemplo trabajaremos con Magna Colombia Bogota EPSG:3116.
Shape de la cuenca
El shape de la cuenca media del Río Sinú junto con las estaciones que serán interpoladas es el que se presenta a continuación:
Alt text
Generación de la malla de puntos
Los centroides de la malla los puedes descargar del drive que he compartido. Sin embargo, sí no sabes generar la malla de puntos, te explicaré de una manera sencilla como hacerlo desde ArcGis.
La malla de puntos se puede generar con la herramienta Fishnet, la cual la puedes encontrar con el botón “Search”:
Alt text
Das click en la primera herramienta que aparece y te debe aparecer la siguiente ventana:
Alt text
Ahora solo hay que diligenciar algunos campos:
-En el espacio Output Feature Class se debe especificar la ruta donde se guardará la malla, para ello se busca la carpeta y se coloca el nombre de tu shape y le das en el botón save.
-En Template Extent (optional) debes dar click y te sale una lista desplegable, luego debes seleccionar el shape de la cuenca.
-Por último debes especificar el tamaño de la celda, para este caso utilizaremos un tamaño de 1500 m x 1500 m, valor que diligenciaremos en el espacio Cell Size Width y Cell Size Height.
Diligenciados los campos te debe quedar algo similar a lo parecido en esta ventana:
Alt text
Con esto ya habremos diligenciado todo lo que se requiere para correr la herramienta Fishnet y daremos click en OK.
Al finalizar te aparecerá una malla cuadriculada con sus respectivos centroides, tal como se muestra a continuación:
Alt text
Ahora solo falta acotar la red de puntos a la cuenca de estudio:
-Para ello puedes hacer uso de la herramienta Select By Location.
-En Target layer(s) seleccionas el shape de los centroides de la malla.
-En Source Layer seleccionas el shape de la cuenca.
Te debe quedar algo parecido a lo siguiente:
Alt text
Das clcick en Apply y listo, te debe aparecer la selección algo parecida a esta:
Alt text
Para exportar los datos seleccionados, das click derecho en el shape de los centroides de la malla y buscas Data y das click en Export Data:
Alt text
Al dar click te aparece la siguiente ventana donde debes asegurarte que en la pestaña deplegable Export diga “Slected features”, buscas la ruta donde quieres guardar y le das OK:
Alt text
Una vez ejecutada la herramienta te debe aparecer una malla de puntos de la siguiente forma:
Alt text
Y listo, ya creaste tu primera malla de puntos que será utilizada como dato de entrada para la interpolación en R.
Datos de precipitacion de las estaciones
Por último requieres de una tabla con los datos de precipitación que puede estar hecha en excel.
El archivo puede o no estar delimitado por comas, para este ejemplo trabajaremos con extensión csv.
Alt text
Haremos interpolación IDW para octubre de 2010.
Código en R
Necesitaremos de las siguientes librerías las cuales deben estar previamente instaladas:
rm(list = ls())
library(xts) # Series Temporales## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(lubridate) # Utilidades para manejo de fechas##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(terra) # Información Geográfica## terra 1.5.21##
## Attaching package: 'terra'## The following object is masked from 'package:zoo':
##
## time<-library(raster)## Loading required package: splibrary(phylin) # IDW
library(rgdal)## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-29, (SVN revision 1165M)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/Usuario/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Usuario/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/Usuario/Documents/R/win-library/4.1/rgdal/proj##
## Attaching package: 'rgdal'## The following object is masked from 'package:terra':
##
## projectSí no sabes instalar una librería en R, lo único que necesitas es darle click en Package en Rstudio y darle click en install:
Alt text
Te aparecerá una ventana como la siguiente:
Alt text
Buscas la librería que necesitas, por ejemplo, para hacer la interpolación mediante el método IDW necesitaremos la librería phylin.
Alt text
Le das click en install y listo!
Mediante la función vect podemos leer el archivo con extensión .shp, aqu[i debes cambiar la ruta que se encuentra entre comillas y adaptarla a la ruta donde se encuentran tus archivos shapes:
estaciones <- vect("C:/Users/Usuario/Documents/4. Parcial/SIG/Est_PP_Totales_3116.shp") # Lee el shp de estacionesPodemos ver parte de la tabla de atributos con el siquiente comando indicando las columnas que queremos visualizar, en este caso tenemos desde la columna 2 a la columna 5 [,2:6]:
as.data.frame(estaciones[,2:6]) # Muestra un extrato de la tabla de atributos## FID_ OBJECTID_2 OBJECTID CÓDIGO NOMBRE
## 1 57 58 646 25010110 ACACIAS LAS HACIENDA [25010110]
## 2 50 51 541 13070170 BUENOS AIRES 1 [13070170]
## 3 34 35 514 13060020 BUENOS AIRES [13060020]
## 4 48 49 537 13070130 CALIFORNIA [13070130]
## 5 47 48 536 13070120 CALLEMAR [13070120]
## 6 3 4 475 13015030 CAMPO BELLO [13015030]
## 7 35 36 515 13060030 CARAMELO [13060030]
## 8 46 47 535 13070110 CARRIZAL [13070110]
## 9 58 59 647 25015010 CENTRO ALEGRE [25015010]
## 10 10 11 484 13025030 DESPENSA LA [13025030]
## 11 30 31 510 13055030 GALAN [13055030]
## 12 51 52 547 13070230 HORIZONTE [13070230]
## 13 19 20 498 13045010 JARAGUAY [13045010]
## 14 52 53 550 13070260 LAMAS 3 [13070260]
## 15 26 27 506 13050030 LOMA VERDE [13050030]
## 16 37 38 517 13065020 MARACAYO [13065020]
## 17 49 50 539 13070150 MOCARI [13070150]
## 18 8 9 568 13075060 MONTERIA [13075060]
## 19 6 7 557 13070340 PALMA DE VINO [13070340]
## 20 18 19 497 13040030 PEZVAL [13040030]
## 21 55 56 643 25010080 PICA PICA [25010080]
## 22 63 64 712 25021550 PLANETA RICA [25021550]
## 23 64 65 724 25025190 PLANETA RICA [25025190]
## 24 12 13 490 13035010 PUERTO NUEVO [13035010]
## 25 17 18 496 13040010 QUIMARI [13040010]
## 26 4 5 476 13015040 REPRESA URRA [13015040]
## 27 54 55 669 25020600 SAJONIA HACIENDA - AUT [25020600]
## 28 33 34 513 13060010 SAN ANTERITO [13060010]
## 29 56 57 645 25010100 SAN FRANCISCO [25010100]
## 30 25 26 505 13050020 SANTA CRUZ HACIENDA [13050020]
## 31 53 54 553 13070290 TAMPA [13070290]
## 32 11 12 489 13030010 TIERRALTA [13030010]
## 33 36 37 516 13065010 TORPEZA LA [13065010]
## 34 61 62 661 25020430 TREMENTINO - AUT [25020430]
## 35 27 28 507 13050040 TRES PIEDRAS [13050040]
## 36 1 2 472 13010020 TUCURA [13010020]
## 37 2 3 474 13015020 TUCURA [13015020]
## 38 7 8 567 13075050 UNIV DE CORDOBA [13075050]
## 39 28 29 508 13055010 VORAGINE LA [13055010]Como te daras cuenta se pueden observar los códigos y nombres de las estaciones, que es lo mismo que tu puedes ver en la tabla de atributos en ArcMap. En este caso interpolaremos 39 estaciones de precipitación.
Mediante la función plot e indicando la columna del campo que tiene los nombres de las estaciones (en este caso 6) podemos visualizar la distribución de las estaciones desde R.
plot(estaciones, 6) # Plotea la ubicación de las estaciones
Necesitamos leer los datos de precipitación de las estaciones para ello utilizaremos la función read.csv, colocamos entre parentesis la ruta del archivo, su nombre y no olvides colocar como extensión .csv.
Utilizamos como separador (sep) el punto y coma ‘;’ esto depende de la configuración de tu computador, debes revisar si la separación es en punto y coma ‘;’ o solo en coma ‘,’.
Datos_PP=read.csv('C:/Users/Usuario/Documents/4. Parcial/Datos_PP_10_2010.csv',sep = ';')Ahora prepararemos algunas variables que se requeriran para el ciclo FOR mediante el cual haremos las interpolaciones de precipitación para 31 días del mes de octubre de 2010.
Asignaremos Datos_PP a otra variable denominada Datos_PP_Interpolados, contaremos las columnas del nuevo datframe mediante la función length y contaremos el número de filas mediante la función nrow.
Datos_PP_Interpolados<-Datos_PP
columnas<-length(Datos_PP_Interpolados)
filas=nrow(Datos_PP_Interpolados)Posteriormente, haremos un ciclo FOR que hara un reccorido desde 1 hasta el número de filas totales, que para este caso, son 31 días.
Para ejecutar el siguiente código, debes cambiar la ruta en la función vect que asigna a la malla nombrada como grid y en path también debes cambiar la ruta a la carpeta propia de tu computador.
Mas ábajo explicaré detalladamente línea por línea, por ahora ejecutaremos el ciclo FOR:
for (i in 1:filas) {
datos.sinteticos <- as.numeric(as.vector(Datos_PP_Interpolados[i,2:40])) # Datos sintéticos
stn <- geom(estaciones)[,c("x", "y")] # Coordenadas de las estaciones
stn <- as.data.frame(stn)
grid <- vect("C:/Users/Usuario/Documents/4. Parcial/SIG/Malla_Puntos_3116_1500.shp") # shp de la red de puntos
grid <- geom(grid)[,c("x", "y")] # coordenadas de los puntos
grid <- as.data.frame(grid)
datos.interpolados <- idw(datos.sinteticos, stn, grid, p = 1) # Aplicación IDW
grid <- data.frame(grid, valores = datos.interpolados[,1])
coordinates(grid) <- ~x+y # Se especifica que 'x' y 'y' son las coordenadas
gridded(grid) <- T # Se indica que se trata de una red de puntos
grid <- raster(grid, "valores") # Se crea el raster
projection(grid) <- crs("+init=epsg:3116") # Se le asigna un sistemas de coordenadas
path='C:/Users/Usuario/Documents/4. Parcial/Raster_Precipitacion/'
numeracion=as.character(i)
Nombre="Campo_PP"
Extension=".tif"
writeRaster(grid, paste(path,Nombre,numeracion,Extension), overwrite = T) # Se escribe el raster
}Si seguiste correctamente el paso a paso en la carpeta de destino tendrás los siguientes archivos raster que puedes visualizar dando doble click en el archivo o también puedes visualizarlos desde ArcMap.
Alt text
Explicación linea a línea el ciclo FOR
Ahora si, explicaré el paso a paso del ciclo FOR para la interpolación del día 1 de octubre de 2010.
La siguiente línea de código lee fila por fila los datos de precipitación de cada estación, por eso i es el contador que aumenta a medida que se avanza la iteración.
Se utiliza la función as.vector para la extracción de los datos y se anida con la función as.numeric para solamente extraer los valores numéricos, los cuales son asignados a la variable datos.sinteticos:
datos.sinteticos <- as.numeric(as.vector(Datos_PP_Interpolados[1,2:40])) # Datos sintéticos La función geom llama al shape de las estaciones de precipitación y extrae los valores de las coordenadas asignandolos a una variable llamada stn:
stn <- geom(estaciones)[,c("x", "y")] # Coordenadas de las estacionesEs necesario convertir la variable stn como dataframe puesto que se requiere que la variable esté en ése tipo formato más adelante:
stn <- as.data.frame(stn)Se carga el shape de la red de puntos mediante la función vect:
grid <- vect("C:/Users/Usuario/Documents/4. Parcial/SIG/Malla_Puntos_3116_1500.shp") Al igual que con la variable stn se extraen las coordenadas de los centroides de la malla y se convierten a tipo dataframe:
grid <- geom(grid)[,c("x", "y")] # coordenadas de los puntos
grid <- as.data.frame(grid)Posteriormente, aplicamos el método de interpolación IDW mediante la función idw que requiere como datos de entrada, los datos que serán interpolados que se encuentran en la variable datos.sinteticos, las coordenadas de la variable datos.sinteticos que se encuentra en la variable stn y los centroides de la malla que serán utilizados en la interpolación:
datos.interpolados <- idw(datos.sinteticos, stn, grid, p = 1) Para este caso se asignó como exponente un número equivalente igual a uno (1) que es lo que usualmente se recomienda para no suavizar los datos en la interpolación.
Seguidamente, se asignan todos los valores interpolados a la malla:
grid <- data.frame(grid, valores = datos.interpolados[,1]) Se asignan las coordenadas a la malla:
coordinates(grid) <- ~x+y # Se especifica que 'x' y 'y' son las coordenadas
gridded(grid) <- T # Se indica que se trata de una red de puntos
grid <- raster(grid, "valores") # Se crea el rasterSe especifica el sistema de coordenadas en el cual estamos trabajando, para este caso es Magna Colombia Bogotá donde su epsg es 3116:
projection(grid) <- crs("+init=epsg:3116") Los últimos comandos se utilizan para indicar la ruta donde se guardará el raster de precipitación. Para ello se crea una variable tipo caracter denominada path que tiene la ruta de la carpeta:
path='C:/Users/Usuario/Documents/4. Parcial/Raster_Precipitacion/'Como i es un iterador numérico se convierte a un caracter y se le asigna a una nueva variable denominada numeración.
i=1
numeracion=as.character(i)Se le da el nombre al archivo raster que se generara y se asigna una variable denominada extensión la cual indica la extensión del archivo raster en formato tif.
Nombre="Campo_PP"
Extension=".tif"Por último, mediante la función writeRaster se guarda la malla generada en formato raster en la carpeta especificada y se utiliza la función paste para concatenar los caracteres y se abilita un proceso de sobreescritura.
writeRaster(grid, paste(path,Nombre,numeracion,Extension), overwrite = T) El overwrite también puede ser False (F) si no se requieren sobreescribir los datos.
Por último, si quieres visualizar tú raster interpolado de precipitación, puedes aplicar el siguiente comando:
plot(grid)
Y listo con esto hemos acabado el tutorial, ya tienes listos tus archivos raster de precipitación en tú carpeta!