Anexo 2: Interpolación poligonos de Thiessen para datos de precipitación en Valle del Cauca

1. Introducción

Este cuaderno se encuentra el anexo número dos del informe final correspondiente a la asignatura geomática básica.En este documento se pretende ilustrar el método de interpolación con polígonos de Thiessen para los datos de precipitación en el Departameno de Vichada

3. Preprocesamiento de datos

Los puntos correspondientes a los datos de precipitación fueron almacenados en un archivo geojson, por ende, es necesario leerlos.

(precip.points <- geojsonio::geojson_read("C:/Users/LUISA CARRION/Documents/chirps/ppoints.geojson11.geojson", what="sp"))

Registered S3 method overwritten by 'geojsonsf':
  method        from   
  print.geojson geojson

class       : SpatialPointsDataFrame 
features    : 677 
extent      : -77.475, -75.725, 3.124999, 5.024999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :           lluvia 
min values  : 17.5390129089355 
max values  : 182.807846069336 

Ahora, se lee el archivo shapefile correspondiente del departamento.

(valle <- shapefile("C:/Users/LUISA CARRION/Downloads/76_VALLE_DEL_CAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 42 
extent      : -77.54977, -75.70724, 3.091239, 5.047394  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                    MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO,      DPTO_CNMBR,     Shape_Leng,       Shape_Area 
min values  :         76,      76001,    ALCALÁ,                          1536,   41.86090736,      2017, VALLE DEL CAUCA, 0.453826056161, 0.00340901158098 
max values  :         76,      76895,     ZARZAL, Ordenanza 9 de Diciembre 1954, 6292.50083741,      2017, VALLE DEL CAUCA,  6.59527269127,   0.510778519244 

Se convierte a un objeto sf empleando la función st_as_sf

valle_sf <-  sf::st_as_sf(valle)
(border_sf <- valle_sf %>%summarise(area = sum(MPIO_NAREA)))

Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
Geodetic CRS:  WGS 84
     area                       geometry
1 20665.6 MULTIPOLYGON (((-77.2346 4....

Se convierte el objeto sf a un Dataframe Espacial.

 (border <- as(border_sf, 'Spatial')) 

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -77.54977, -75.70724, 3.091239, 5.047394  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :           area 
value       : 20665.59885853 

Al objeto creado y denominado VAlle_sf se le agregan las nuevas columnas “MUNIC” y “CODIGO”, que contienen la misma información que “MPIO_CNMBR” y “MPIO_CCDGO” respectivamente.

(valle.sf <- st_as_sf(valle) %>% mutate(MUNIC = MPIO_CNMBR, CODIGO = MPIO_CCDGO) %>% select(MUNIC, CODIGO))

Simple feature collection with 42 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
Geodetic CRS:  WGS 84
First 10 features:
         MUNIC CODIGO                       geometry
0         CALI  76001 MULTIPOLYGON (((-76.59175 3...
1 ANDALUCÃ\u008dA  76036 MULTIPOLYGON (((-76.22406 4...
2 ANSERMANUEVO  76041 MULTIPOLYGON (((-76.01558 4...
3      ARGELIA  76054 MULTIPOLYGON (((-76.14316 4...
4         BUGA  76111 MULTIPOLYGON (((-76.31608 3...
5 BUGALAGRANDE  76113 MULTIPOLYGON (((-76.15131 4...
6   CAICEDONIA  76122 MULTIPOLYGON (((-75.8539 4....
7       CALIMA  76126 MULTIPOLYGON (((-76.51747 4...
8   CANDELARIA  76130 MULTIPOLYGON (((-76.30455 3...
9      CARTAGO  76147 MULTIPOLYGON (((-75.94518 4...

Los puntos de datos de precipitación se convierten en un objeto sf

p.sf <- st_as_sf(precip.points)

Empleando la función “st_intersect” se realiza una interseccióon de los poligonos con los puntos.

(precip.sf = st_intersection(valle.sf, p.sf))

Warning: attribute variables are assumed to be spatially constant throughout all geometries

Simple feature collection with 677 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -77.475 ymin: 3.124999 xmax: -75.725 ymax: 5.024999
Geodetic CRS:  WGS 84
First 10 features:
            MUNIC CODIGO   lluvia                 geometry
10   EL Ã\u0081GUILA  76243 59.29984 POINT (-76.075 5.024999)
10.1 EL Ã\u0081GUILA  76243 65.32109 POINT (-76.075 4.974999)
10.2 EL Ã\u0081GUILA  76243 68.39664 POINT (-76.075 4.924999)
10.3 EL Ã\u0081GUILA  76243 65.94273 POINT (-76.025 4.924999)
10.4 EL Ã\u0081GUILA  76243 64.68098 POINT (-76.125 4.874999)
10.5 EL Ã\u0081GUILA  76243 66.65178 POINT (-76.075 4.874999)
10.6 EL Ã\u0081GUILA  76243 64.36810 POINT (-76.025 4.874999)
11       EL CAIRO  76246 71.65340 POINT (-76.175 4.824999)
2    ANSERMANUEVO  76041 58.66259 POINT (-76.125 4.824999)
2.1  ANSERMANUEVO  76041 63.25397 POINT (-76.075 4.824999)

Ahora, se hace una reproyección en ambos objetos.

p.sf.magna <- st_transform(precip.sf, crs=3116)
valle.sf.magna <- st_transform(valle.sf, crs=3116)

Se realiza una conversión en el objeto “p.sf.magma”

(precip2 <- as(p.sf.magna, 'Spatial'))

class       : SpatialPointsDataFrame 
features    : 677 
extent      : 622166.5, 817183.3, 837632.5, 1047756  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 3
names       :   MUNIC, CODIGO,           lluvia 
min values  : ALCALÁ,  76001, 17.5390129089355 
max values  :  ZARZAL,  76895, 182.807846069336 

Se guarda el conjunto de datos obtenido en la carpeta chirps con el fin de tener una copia que pueda ser usada después.

shapefile(precip2, filename='chirps/precip2.shp', overwrite=TRUE)

Ahora, se redondearán los valores de la variable “lluvia” y dichos resultados, se almacenaran en la columna “precipitación”

precip2$precipitacion <- round(precip2$lluvia, 1)
precip2

class       : SpatialPointsDataFrame 
features    : 677 
extent      : 622166.5, 817183.3, 837632.5, 1047756  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       :   MUNIC, CODIGO,           lluvia, precipitacion 
min values  : ALCALÁ,  76001, 17.5390129089355,          17.5 
max values  :  ZARZAL,  76895, 182.807846069336,         182.8 

Posteriormente, se realiza una conversión al objeto sf (“valle.sf.magma”) ,

(valle2 <- as(valle.sf.magna, 'Spatial'))

class       : SpatialPolygonsDataFrame 
features    : 42 
extent      : 613842.8, 819150.3, 833918.1, 1050235  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :   MUNIC, CODIGO 
min values  : ALCALÁ,  76001 
max values  :  ZARZAL,  76895 

shapefile(valle2, filename='chirps/valle2.shp', overwrite=TRUE)

Con el fin de verificar que las dos extensiones coincidan, se reemplaza la extensión del limite de los puntos con la del objeto “valle 2”

precip2@bbox <- valle2@bbox

Finalmente, se grafican los datos obtenidos luego del preprocesamiento

tm_shape(valle2) + tm_polygons() +tm_shape(precip2) +tm_dots(col="precipitacion", palette = "Dark2", midpoint = 42, title="Muestras de precipitacion \n(mm)", size=0.08) + tm_legend(legend.outside=TRUE)

Warning in sp::proj4string(obj) :
  CRS object has comment, which is lost in output

4. Interpolación con polígonos de Thiessen

Se representa el patron de puntos bidimensionales como formato ppp para que se pueda usar la función dirichlet(). Esta función es la que permite crear la superficie teselada de los polígonos de Thiessen.

library(spatstat)
th  <-  as(dirichlet(as.ppp(precip2)), "SpatialPolygons")
crs(th) <- crs(precip2)
crs(valle2) <- crs(precip2)

Se emplea la función over() para unir los atributos del punto a la superficie mediante una unión espacial. Tambien es necesario recortar la superficie al tamaño del departamento.

th.z     <- over(th, precip2, fn=mean)
th.spdf  <-  SpatialPolygonsDataFrame(th, th.z)
th.clp   <- raster::intersect(valle2,th.spdf)

Se grafica la información obtenida.

tm_shape(th.clp) + tm_polygons(col="precipitacion", palette="Spectral", midpoint=43.0, title="Polígonos de Thiessen \n para precipitación \n(mm)") + tm_legend(legend.outside=TRUE,title.size=1.2, text.size= 0.8)

Warning in sp::proj4string(obj) :
  CRS object has comment, which is lost in output

5. Otra forma de interpolación con poígonos de Thiessen

Se toma aleatoramente la mitad de los datos del objeto “precip2” con la función “sample()”.

train_index <- sample(1:nrow(precip2), 0.7 * nrow(precip2))
test_index <- setdiff(1:nrow(precip2), train_index)
ptos_train <- precip2[train_index, ]
ptos_test  <- precip2[test_index,]

ptrain <- spTransform(ptos_train, crs(precip.mask))
ptest <- spTransform(ptos_test, crs(precip.mask))

Se crea un mapa que permite de la visualización de los 2 conjuntos de datos.

library(htmltools)
library(leaflet.extras)

lplot <- leaflet(data = precip2) %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addRasterImage(precip.mask, colors = pal, opacity = 0.6) %>%
  addCircleMarkers(data = ptrain, 
                   radius = 1,
                   fillOpacity = .7,
                   stroke = FALSE,
                   popup = ~htmlEscape(precipitacion),
                   color = pal(ptos_train$precipitacion), 
                   clusterOptions = markerClusterOptions(),
                   group = "Training") %>% 
  addCircleMarkers(data = ptest, # second group
                   radius = 10,
                   fillOpacity = .7,
                   stroke = FALSE,
                   popup = ~htmlEscape(precipitacion),
                   color = pal(ptos_test$precipitacion), 
                   clusterOptions = markerClusterOptions(),
                   group = "Test") %>% 
  addLegend(position = "bottomright",
            values = ~precipitacion,
            opacity = .7,
            pal = pal, 
            title = "Precipitacion") %>% 
  leaflet::addLayersControl(overlayGroups = c("Training", "Test"),
                   options = layersControlOptions(collapsed = FALSE)) %>% 
  addResetMapButton()

lplot

Se crea la superficie teselada.

thiessen  <-  as(dirichlet(as.ppp(ptos_train)), "SpatialPolygons")

Se añade la información de proyección.

crs(thiessen) <- crs(ptos_train)
crs(valle2) <- crs(ptos_train)

Se unen los atributos

thiessen.z     <- over(thiessen, ptos_train, fn=mean)
thiessen.spdf  <-  SpatialPolygonsDataFrame(thiessen, thiessen.z)

Se recorta la superficie al tamaño del departamento.

thiessen.clp   <- raster::intersect(valle2,thiessen.spdf)

Se grafican los polígonos de Thiessen del departamento del Valle del Cauca.

tm_shape(thiessen.clp) + tm_polygons(col="precipitacion", palette="RdYlGn" , midpoint=43.0, title="Poligonos de Thiessen \npara precipitación \n(mm)") + tm_legend(legend.outside=TRUE, title.size=1.2, text.size= 0.8)

Warning in sp::proj4string(obj) :
  CRS object has comment, which is lost in output