Interpolación poligonos de Thiessen para datos de precipitación en Vichada

1. Introducción

Este cuaderno , correspone al anexo número dos del informe final correspondiente a la asignatura geomática básica .En el se ilustra 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 en el departamento de Vichada , fueron obtenidos y almacenados en un archivo geojson.Ahora , usando la función geojson_read() , se leen estos puntos

(precip.points <- geojsonio::geojson_read("./chirps/ppoints.geojson", what="sp"))

class       : SpatialPointsDataFrame 
features    : 3248 
extent      : -71.075, -67.475, 2.774999, 6.274999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       :            lluvia 
min values  : 0.541795015335083 
max values  :  85.8214721679688 

There were 24 warnings (use warnings() to see them)

Se lee de nuevo el archivo shapefile correspondiente al área de estudio

(Vichada <- shapefile("./ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 4 
extent      : -71.07793, -67.4098, 2.737109, 6.324317  (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  :         99,      99001,      CUMARIBO,       Decreto 1594 de Ago 5  de 1974, 3898.56891769,      2017,    VICHADA, 3.29670807195, 0.316732435778 
max values  :         99,      99773, SANTA ROSALÍA, Ordenanza 66 de Noviembre 22 de 1996, 65599.7022767,      2017,    VICHADA,  18.794382661,   5.3085802966 

Posteriormente, se convierte a un objeto sf empleando la función st_as_sf y se disuelven los límites internos como se ilustra en el siguiente código

Vichada_sf <-  sf::st_as_sf(Vichada)
(border_sf <- Vichada_sf %>%summarise(area = sum(MPIO_NAREA)))

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -71.07793 ymin: 2.737109 xmax: -67.4098 ymax: 6.324317
CRS:            +proj=longlat +datum=WGS84 +no_defs 
      area                       geometry
1 100065.9 POLYGON ((-67.814 5.342584,...

Se convierte el objeto sf a un Dataframe Espacial . En este caso, se refiere a un polígono que contiene atributos

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

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -71.07793, -67.4098, 2.737109, 6.324317  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       :            area 
value       : 100065.86486549 

Al objeto Vichadasf , se le agregan las nuevas columnas “MUNIC” y “CODIGO” , que tienen el mismo contenido que “MPIO_CNMBR” y “MPIO_CCDGO” y debido a que estas nuevas columnas se agregan al final del conjunto de datos , se usa la función select() para que solo se vean estas variables en la salida del código

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

Simple feature collection with 4 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -71.07793 ymin: 2.737109 xmax: -67.4098 ymax: 6.324317
CRS:            +proj=longlat +datum=WGS84 +no_defs 
           MUNIC CODIGO                       geometry
1  SANTA ROSALÍA  99624 POLYGON ((-70.65378 5.37297...
2 PUERTO CARREÑO  99001 POLYGON ((-67.80972 6.32431...
3   LA PRIMAVERA  99524 POLYGON ((-69.03359 6.21869...
4       CUMARIBO  99773 POLYGON ((-68.47074 5.55046...

Los puntos de datos de precipitación , también son convertidos a un objeto sf

p.sf <- st_as_sf(precip.points)

Usando la función st_intersect de la libreria sf, se realizará una intersección de los polígonos con los puntos, manteniendo la información de los 2

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

although coordinates are longitude/latitude, st_intersection assumes that they are planar

Simple feature collection with 3248 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -71.075 ymin: 2.774999 xmax: -67.475 ymax: 6.274999
CRS:            +proj=longlat +datum=WGS84 +no_defs 
First 10 features:
             MUNIC CODIGO    lluvia                 geometry
2   PUERTO CARREÑO  99001 0.8161542 POINT (-67.875 6.274999)
2.1 PUERTO CARREÑO  99001 0.8258479 POINT (-67.825 6.274999)
2.2 PUERTO CARREÑO  99001 0.8074083 POINT (-67.775 6.274999)
2.3 PUERTO CARREÑO  99001 0.8239177 POINT (-67.725 6.274999)
2.4 PUERTO CARREÑO  99001 0.7755996 POINT (-67.675 6.274999)
2.5 PUERTO CARREÑO  99001 0.7823746 POINT (-67.625 6.274999)
2.6 PUERTO CARREÑO  99001 3.5164738 POINT (-68.125 6.224999)
2.7 PUERTO CARREÑO  99001 3.2473364 POINT (-68.075 6.224999)
2.8 PUERTO CARREÑO  99001 0.7987307 POINT (-67.925 6.224999)
2.9 PUERTO CARREÑO  99001 0.7961798 POINT (-67.875 6.224999)

Se realiza una tarea de reproyección en los 2 objetos

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

Luego, se hace una conversión del objeto “p.sf.magma”

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

class       : SpatialPointsDataFrame 
features    : 3248 
extent      : 1333188, 1732620, 799126.6, 1190067  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +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  :      CUMARIBO,  99001, 0.541795015335083 
max values  : SANTA ROSALÍA,  99773,  85.8214721679688 

Se puede guardar el conjunto de datos en la carpeta chirps del directorio de trabajo para que pueda ser usado de forma posterior

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

Ahora, se aproximan los valores que toma la variable “lluvia” a un decimal.A la columna que tiene estos valores redondeados, se le asignará el nombre “precipitación”

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

class       : SpatialPointsDataFrame 
features    : 3248 
extent      : 1333188, 1732620, 799126.6, 1190067  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +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  :      CUMARIBO,  99001, 0.541795015335083,           0.5 
max values  : SANTA ROSALÍA,  99773,  85.8214721679688,          85.8 

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

(Vichada2 <- as(Vichada.sf.magna, 'Spatial'))

class       : SpatialPolygonsDataFrame 
features    : 4 
extent      : 1332850, 1739856, 794950.2, 1195302  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +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  :      CUMARIBO,  99001 
max values  : SANTA ROSALÍA,  99773 

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

Para asegurarse de que las dos extensiones coincidan se reemplaza la extensión del límite de los puntos con la del nuevo objeto “Vichada 2”

precip2@bbox <- Vichada2@bbox

Se plotean los datos obtenidos luego del preprocesamiento

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

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() de la libreria spatstat .Esta función es la que permite crear la superficie teselada de los polígonos de Thiessen , la superficie generada no tiene información de proyección por tanto, también se añade de forma manual

th  <-  as(dirichlet(as.ppp(precip2)), "SpatialPolygons")
crs(th) <- crs(precip2)
crs(Vichada2) <- crs(precip2)

Como , no almacena información de atributos de la capa de puntos. se usa la función over() (del paquete sp) para unir los atributos del punto a la superficie mediante una unión espacial. Posteriormente , se recorta la superficie al tamaño del departamento

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

Ahora, se plotea la superficie interpolada

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

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

Con la función sample(),se toma una muestra aleatoria de la mitad los datos del objeto “precip2” .Para obtener los datos de “precip2” que no están en la variable creada , se hace uso de la función setdiff()

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

Se crea un mapa que permite de la visualización de los 2 conjuntos de datos con ayuda de las librerias leaflet, htmltools y leaflet.extras

# Se agrega la capa ráster de datos de precipitación 
lplot <- leaflet(data = precip2) %>% addProviderTiles("CartoDB.Positron") %>% addRasterImage(precip.mask, colors = pal, opacity = 0.6, project=FALSE) %>% 
# Se agrega el primer conjunto de datos
addCircleMarkers(data = ptrain, radius = 1, fillOpacity = .7, stroke = FALSE, popup = ~htmlEscape(precipitacion), color = pal(ptos_train$precipitacion), clusterOptions = markerClusterOptions(),group = "Training") %>% 
# Se agrega el segundo conjunto de datos  
addCircleMarkers(data = ptest, radius = 10, fillOpacity = .7, stroke = FALSE, popup = ~htmlEscape(precipitacion), color = pal(ptos_test$precipitacion), clusterOptions = markerClusterOptions(), group = "Test") %>% 
# Se ajustan las opciones de la leyenda y las opciones para activar o desactivar las capas  
addLegend(position = "bottomright",values = ~precipitacion, opacity = .7, pal = pal, title = "Precipitacion") %>% leaflet::addLayersControl(overlayGroups = c("Training", "Test"), options = layersControlOptions(collapsed = FALSE)) %>% addResetMapButton()

lplot

Se reproyectan los datos seleccionados

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

Se crea la superficie teselada , luego de convertir los datos a formato ppp

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

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

crs(thiessen) <- crs(ptos_train)
crs(Vichada2) <- crs(ptos_train)

Se unen los atributos a la superficie teselada

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

Se recortan la superficie al tamaño de la zona de interés

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

Se plotean los polígonos de thiessen obtenidos para el departamento de Vichada

tm_shape(thiessen.clp) + tm_polygons(col="precipitacion", palette="RdBu" , 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)