OBJETIVO

La puesta en práctica de este script tiene como finalidad la estimación de la humedad del suelo a partir de datos obtenidos por teledetección. El área de estudio se localiza en el departamento de Quindío en los municipios de “ALCALÁ-ARMENIA-BUENAVISTA-CAICEDONIA-CALARCÁ-CIRCASIA-CÓRDOBA-FILANDIA-GÉNOVA- LA TEBAIDA-MONTENEGRO-PIJAO-QUIMBAYA-SALENTO-ULLOA”

DATOS

Se utilizan dos fuentes diferentes para obtener variables explicativas:

SENTINEL 1-A

  1. un conjunto de datos de imágenes de radar de la misión Sentinel-1A con un sensor activo de radar de apertura sintética (SAR), proporciona imágenes en banda C tanto en polarización singular VV como en polarización dual VH, no se ve afectado por la cobertura de nubes o la falta de iluminación. Las imágenes son de acceso libre a cualquier usuario, las imágenes s de adquisición son en modo interferométrico de amplia franja (IW), permite la combinación de una anchura de barrido de gran tamaño (250 km) con una resolución geométrica moderada de 5 m x 20 m. La fecha de adquisición de los datos es de Mayo 12 de 2015 con una resolucion de 20 metros.

PORTAL ESA

PORTAL ESA COPERNICUS

PORTAL ESA “COPERNICUS”

SENTINEL 2-A

  1. imágenes multiespectrales SENTINEL 2A y su sensor MSI (Multi Spectral Imager) de resolución espacial alta-media, caracterizados por 290 kilómetros ancho de franja y una alta capacidad de revisita (5 días con dos satélites), potencialmente adecuado para el mapeo regional de coberturas de la tierra. Este sensor ofrece un diseño polivalente de 13 bandas espectrales que atraviesan desde la región visible del espectro e infrarrojo cercano hasta el infrarrojo de onda corta, Sentinel 2 con su sensor MSI cuenta con cuatro bandas (2, 3, 4 y 8) con una resolución espacial de 10 m, seis bandas (5, 6, 7, 8a, 11 y 12) a 20m y las últimas tres bandas (1, 9 y 10) a 60 m, entre estas trece bandas hay tres nuevas bandas en la región del rojo posicionadas a 705, 740 y 783 nm. De estas imágenes se calcularon cuatro índices que fueron dos NDVI implementando dos infrarrojos NIR (B5 - B8) y dos índices de diferencia normalizada de agua NDWI implementando dos infrarrojos (B5 - B8) y el SWIR (B11) como se resume en la tabla.La fecha de adquisición de los datos es de Ocubtre 01 de 2015 con una resolucion de 20 metros.
BANDAS SENTINEL 2A

BANDAS SENTINEL 2A

DATOS DE HUMEDAD DEL SUELO SMAP

La propiedad del suelo objetivo de estudio es la humedad del suelo superficial. Estos datos provienen del geoportal Monitor de inundaciones y sequía en América Latina y el Caribe (Latin American and Caribbean Flood and Drought Monitor “LAFDM”) donde se encuentran datos de humedad del suelo de la misión SMAP (Soil Moisture Active Passive); es una misión pasiva activa de humedad del suelo que utiliza un radar de banda L y un radiómetro de banda L para mediciones concurrentes y coincidentes integradas como un solo sistema de observación. Esta combinación aprovecha las ventajas relativas de la detección remota de microondas tanto activa (radar) como pasiva (radiómetro) para el mapeo de humedad del suelo. La fecha de adquisición de los datos es de Mayo 01 de 2015 escalados a una resolucion de 20 metros implementando un modelo de geosetadistica KRIGING CIRCULAR.

PORTAL LAFDM

BANDAS SENTINEL 2A

BANDAS SENTINEL 2A

ÁREA DE ESTUDIO

El area de studio se comopone de 15 municipios “ALCALÁ-ARMENIA-BUENAVISTA-CAICEDONIA-CALARCÁ-CIRCASIA-CÓRDOBA-FILANDIA-GÉNOVA- LA TEBAIDA-MONTENEGRO-PIJAO-QUIMBAYA-SALENTO-ULLOA” que se pueden ver en la siguiente figura.

#### LIBRERIAS REQUERIDAS
library(raster)
library(rasterVis)
library(tmap)
library(rgdal)
library(randomForest)

Loading required package: sp
Loading required package: lattice
Loading required package: latticeExtra
Loading required package: RColorBrewer

setwd("C:/Users/Usuario/Desktop/PRA/UTIMO TRABAJO PRA/New Folder")
list.files(".",pattern ="shp")
[1] "MUN.shp"     "PUNTOS_MOL.shp"  "PUNTOS_SMAP.shp"  "AE_II.shp"  "quindioII.shp" 

AE<- shapefile("AE_II.shp")
Q<-shapefile("quindioII.shp")

m <- leaflet(Q) %>% 
  addPolygons() %>% 
  addTiles()
  addMiniMap(m, tiles =  providers$Esri.WorldStreetMap,
             toggleDisplay = T, width = 120, height=100)

DEPARTAMENTO DE QINDIO

DEPARTAMENTO DE QUINDIO DEPARTAMENTO DE QUINDIO

ZONA DE ESTUDIO

ZONA DE ESTUDIO

n <- leaflet(AE) %>% 
  addPolygons() %>% 
  addTiles()
  addMiniMap(n, tiles =  providers$Esri.WorldStreetMap,
             toggleDisplay = T, width = 120, height=100)
ZONA DE ESTUDIO

ZONA DE ESTUDIO

DATOS DE HUMEDAD DEL SUELO

Los datos de humedad provienen del portal LAFDM de la misión SMAP con una resolución espacial de 30 km. Estos datos fueron escalados a 20 metros implementando un kriging circular

list.files(".",pattern = "tif")

P<-shapefile("quindio.shp")

HUM_SMAP=raster("smap_SOIL_M.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(HUM_SMAP, margin=F, par.settings=mapTheme, at = do.breaks(c(0,100),10),
                 main="SMAP PORCENTAJE DE HUMEDAD PARA COLOMBIA")
plt + layer(sp.lines(P, col="red", lwd=0.5))
DATOS SMAP

DATOS SMAP

DATOS DE % HUMEDAD ESCALADOS AL DEPARTAMENTO DE QUINDIO Y ZONA DE ESTUDIO

AE_I<-shapefile("ZONA_E.shp")

HUM_QUINDIO=raster("HUM_III.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(HUM_QUINDIO, margin=F, par.settings=mapTheme, at = do.breaks(c(70,100),10),
                 main="PORCENTAJE DE HUMEDAD ESCALADO ZONA DE ESTUDIO")
plt + layer(sp.lines(AE_I, col="black", lwd=0.9))
DATOS SMAP ESCALADOS

DATOS SMAP ESCALADOS

DATOS SENTINEL SENTINEL 2-A

NDVI NIR (B5)

AE_I<-shapefile("ZONA_E.shp")

NDVI_B5=raster("NDVI_B5.tif")

mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
 
plt <- levelplot(NDVI_B5, margin=F, par.settings=mapTheme, at = do.breaks(c(-1,1),30),
+                  main="NDVI_B5 ZONA DE ESTUDIO")
plt + layer(sp.lines(AE_I, col="black", lwd=1))
NDVI_B5

NDVI_B5

NDVI NIR (B8)

AE_I<-shapefile("ZONA_E.shp")

NDVI_B8=raster("NDVI_B8.tif")

mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
 
plt <- levelplot(NDVI_B8, margin=F, par.settings=mapTheme, at = do.breaks(c(-1,1),30),main="NDVI_B8 ZONA DE ESTUDIO")
plt + layer(sp.lines(AE_I, col="black", lwd=1))
NDVI_B8

NDVI_B8

NDWI NIR (B5)

AE_I<-shapefile("ZONA_E.shp")

NDWI_B5=raster("NDWI_B5.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(NDWI_B5, margin=F, par.settings=mapTheme, at = do.breaks(c(-1,1),30),
                 main="NDWI_B5 ZONA DE ESTUDIO")
plt + layer(sp.lines(AE_I, col="black", lwd=0.9))
NDWI_B5

NDWI_B5

NDWI NIR (B8)

AE_I<-shapefile("ZONA_E.shp")

NDWI_B8=raster("NDWI_B8.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(NDWI_B8, margin=F, par.settings=mapTheme, at = do.breaks(c(-1,1),30),
                 main="NDWI_B8 ZONA DE ESTUDIO")
plt + layer(sp.lines(AE_I, col="black", lwd=0.9))
NDWI_B8

NDWI_B8

DATOS SENTINEL SENTINEL 1-A

Las imagen Sentinel-1A fueron calibradas en el ángulo de incidencia Sigma, fue aplicado el proceso de “multilooking” a dos rangos con el fin de obtener pixeles cuadrados, se corrigió el “Speckle” o ruido con el filtro Lee Sigma en una ventana de 5x5 y fue corregida sobre terreno con el modelo digital de elevación SRTM 1sec. Los valores de niveles digitales (DN) de la imagen SAR fueron convertidos a valores de retrodispersión “Backscattering” en la escala de decibelios (db), estos procesos se llevaron a cabo en el software SNAP, el cual es proporcionado por la ESA en la plataforma de aplicaciones Sentinel.

PAGINA ESA DESCARGA SOFTWARE SNAP

BANDA VV

AE_I<-shapefile("ZONA_E.shp")

VV=raster("VV.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(VV, margin=F, par.settings=mapTheme, at = do.breaks(c(-18,24),20),
                 main="IMAGEN DE RADAR POLARIZACION VV ZONA DE ESTUDIO")
plt + layer(sp.lines(AE_I, col="black", lwd=0.9))
IMAGEN DE RADAR POLARIZACION VV

IMAGEN DE RADAR POLARIZACION VV

BANDA VH

AE_I<-shapefile("ZONA_E.shp")

VH=raster("VH.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(VH, margin=F, par.settings=mapTheme, at = do.breaks(c(-23,3),20),
                 main="IMAGEN DE RADAR POLARIZACION VH ZONA DE ESTUDIO")
plt + layer(sp.lines(AE_I, col="black", lwd=0.9))
IMAGEN DE RADAR POLARIZACION VH

IMAGEN DE RADAR POLARIZACION VH

IMPLEMENTACIÓN DEL ALGORITMO DE APRENDIZAJE RANDOMFOREST

Fueron seleccionados 180 puntos aleatoriamente sobre los cuales se extrajo el valor de la variable de interés HUMEDAD DEL SUELO de los datos escalados de SMAP en la zona de estudio como se muestra en la imagen con el fin de entrenar el algoritmo RandomForest.

PUNTOS<-shapefile("PUNTOS_MOL.shp")

HUM=raster("HUM.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(HUM, margin=F, par.settings=mapTheme, at = do.breaks(c(75,82),20),
                 main="180 PUNTOS DE ENTRENAMIENTO")
plt + layer(sp.lines(PUNTOS, col="black", lwd=2))
PUNTOS DE ENTRENAMIEMTO PARA EL ALGORITMO RANDOM FOREST

PUNTOS DE ENTRENAMIEMTO PARA EL ALGORITMO RANDOM FOREST

variables <- list.files(pattern='.tif') #estraer las variables raster
variables
[1] "NDVI_B5.tif" "NDVI_B8.tif" "NDWI_B5.tif" "NDWI_B8.tif" "VH.tif"     
[6] "VV.tif"  


var_raster <- raster(variables[1])#cargar cualquier variable
spplot(var_raster)
VARIABLE 1 NDVI_B5

VARIABLE 1 NDVI_B5

rbrick <- var_raster  #asignar a la variable rbrick
for (i in 2:length(variables)) #leer el resto de variables
{ var_raster <- raster(variables[i])
rbrick <- addLayer(rbrick,var_raster)
}

rbrick

class       : RasterStack 
dimensions  : 842, 900, 757800, 6  (nrow, ncol, ncell, nlayers)
resolution  : 20, 20  (x, y)
extent      : 1135866, 1153866, 980483, 997323  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
names       : NDVI_B5, NDVI_B8, NDWI_B5, NDWI_B8, VH, VV 

names(rbrick)

[1] "NDVI_B5" "NDVI_B8" "NDWI_B5" "NDWI_B8" "VH"  "VV"   

summary(rbrick)

 NDVI_B5     NDVI_B8    NDWI_B5    NDWI_B8        VH         VV
Min.    -0.4573770 -0.04005086 -0.5006797 -0.5590400 -22.56884 -16.859678
1st Qu.  0.1518518  0.52439156 -0.3420253  0.1100038 -14.17802  -8.024954
Median   0.2348050  0.64572763 -0.3106405  0.2087868 -13.46359  -7.000395
3rd Qu.  0.2920257  0.71383770 -0.2649221  0.2846662 -12.77365  -5.910652
Max.     0.7062821  0.83405387  0.3630769  0.7599838   0.00000  17.328157
NA's     0.0000000  0.00000000  0.0000000  0.0000000   0.00000   0.000000

spplot(rbrick)

VARIABLES PREDICTORAS

LEER PUNTOS DE ENTRENAMIENTO
entren<-readOGR(".","PUNTOS_MOL")

OGR data source with driver: ESRI Shapefile 
Source: ".", layer: "PUNTOS_MOL"
with 180 features
It has 2 fields
Integer64 fields read as strings:  CID 
head(entren@data)

 CID  SOIL_M
1   0 80.2631
2   0 79.4950
3   0 77.4009
4   0 79.2618
5   0 76.4223
6   0 77.8940

RASTERIZAR PUNTOS

SM <- rasterize(entren, rbrick, 'SOIL_M')
spplot(SM, scales=list(draw=T))

PUNTOS RASTERIZADOS

CREAR MASCARA SOBRE LAS VARIABLES PREDICTORAS Y VARIABLE DEPENDIENTE

mascara <- mask(rbrick, SM)
entren_RF <- addLayer(mascara, SM)
plot(entren_RF)

MASCARA CON TODAS LAS VARIABLES

RENOMBRAMOS LAS VARIABLES

names(entren_RF)<- c("NDVI_B5", "NDVI_B8", "NDWI_B5", "NDWI_B8", "VH","VV","SM_HUM")

GENERAMOS LA TABLA DE VALORES

tabVAL <- getValues(entren_RF)
tabVAL <- na.omit(tabVAL)
tabVAL <- as.data.frame(tabVAL)
tabVAL$clase <- factor(tabVAL$SM_HUM, levels = c(1:180))

ENTRENAMOS EL ALGORITMO

modelRF <- randomForest(x=tabVAL[ ,c(1:6)], y=tabVAL$SM_HUM ,ntree=1000, importance = TRUE, do.trace=FALSE)

LLAMAMOS AL MODELO PARA VERIFICAR EL PORCENTAJE DE EXPLICACION Y EL CUADRADO MEDIO DE RESIDUALES

modelRF 
Call:
 randomForest(x = tabVAL[, c(1:6)], y = tabVAL$SM_HUM, ntree = 1000,      importance = TRUE, do.trace = FALSE) 
               Type of random forest: regression
                     Number of trees: 1000
No. of variables tried at each split: 2

          Mean of squared residuals: 1.599992
                    % Var explained: 8.75

EVALUAMOS LA IMPORTANCIA DE CADA PREDICTOR

importance(modelRF)

          %IncMSE IncNodePurity
NDVI_B5  9.461432      48.79808
NDVI_B8  8.869579      41.77878
NDWI_B5 17.316276      55.47305
NDWI_B8  6.376657      41.02008
VH       9.969921      48.76350
VV      18.446685      55.15847

GRAFICAMOS LA IMPORTANCIA DE CADA PREDICTOR

varImpPlot(modelRF)

APLICAMOS EL ALGORITMO ENTRENADO A LAS VARIABLES PREDICTORAS

pred_SM <- predict(rbrick, model=modelRF, filename="predicccion_SM.tif", na.rm=TRUE, overwrite=TRUE)
save(modelRF, file="prediction_modelRF.RData")

PLOTEAMOS EL RESULTADO DE LA PREDICCION

AE_I<-shapefile("ZONA_E.shp")

PREDICCION=raster("predicccion_SM.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(PREDICCION, margin=F, par.settings=mapTheme, at = do.breaks(c(75,80),40),
                 main="PREDICCION DE HUMEDAD RF")
plt + layer(sp.lines(AE_I, col="black", lwd=0.9))

GENERAMOS UN RESUMEN EXPLORATORIO DE LA PREDICCIÓN

summary(PREDICCION)

predicccion_SM
Min.          76.34333
1st Qu.       77.82462
Median        78.22325
3rd Qu.       78.53077
Max.          79.72274
NA's           0.00000

COMPARAMOS LOS VALORES DE HUMEDAD “SMAP” CON LOS VALORES PREDICHOS “RANDOMFOREST”

AE_I<-shapefile("ZONA_E.shp")

HUM=raster("HUM.tif")

mapTheme <- rasterTheme(region=brewer.pal(15,"Spectral"))

plt <- levelplot(HUM, margin=F, par.settings=mapTheme, at = do.breaks(c(75,82),20),
                 main="% HUMEDAD SMAP")
plt + layer(sp.lines(AE_I, col="black", lwd=2))
AE_I<-shapefile("ZONA_E.shp")

PREDICCION=raster("predicccion_SM.tif")

mapTheme <- rasterTheme(region=brewer.pal(11,"Spectral"))

plt <- levelplot(PREDICCION, margin=F, par.settings=mapTheme, at = do.breaks(c(75,80),20),
                 main="PREDICCION DE HUMEDAD RF")
plt + layer(sp.lines(AE_I, col="black", lwd=0.9))

GRAFICAMOS EL HISTOGRAMA DE EL RASTER DE HUMEDAD DE ENTRENAMIENTO “DATOS SMAP” Y EL RASTER DE PREDICCION

par(mfrow=c(1,1))

hist(HUM,col="blue", xlab="%HUMEDAD",  main="DATOS SMAP")

hist(PREDICCION,col="green",xlab="%HUMEDAD", main="PREDICCION DE HUMEDAD")

REALIZAMOS UNA DIFERENCIA ENTRE EL VALOR OBSERVADO Y EL PREDICHO

DIF <-(HUM - PREDICCION)

REALIZAMOS UN RESUMEN EXPLORATORIO DE LA DIFERENCIA CON SU RESPECTIVO HISTOGRAMA

summary(DIF)

layer
Min.    -3.7104263
1st Qu. -0.8756142
Median   0.1176758
3rd Qu.  1.0978775
Max.     3.9222412
NA's     0.0000000

hist(DIF,col="orange", xlab="HUMEDAD",  main="RESIDUAL")

GRAFICAMOS EL RASTER DE LA DIFERENCIA

AE_I<-shapefile("ZONA_E.shp")

mapTheme <- rasterTheme(region=brewer.pal(5,"Spectral"))

plt <- levelplot(DIF, margin=F, par.settings=mapTheme, at = do.breaks(c(-4,4),3),
                 main="RASTER DE RESIDUALES")
plt + layer(sp.lines(AE_I, col="black", lwd=0.9))
EN AZUL LAS ZONAS DE SOBREESTIMACION DE LA PREDICCION SUPERIORES A 1% DE HUMEDAD , EN ROJO LAS ZONAS DE SUBESTIMACION INFERIORES A 1% DE HUMEDAD Y EN COLOR CLARO LAS ZONAS CON PREDICCIONES ENTRE -1% Y 1% DE HUEMDAD

EN AZUL LAS ZONAS DE SOBREESTIMACION DE LA PREDICCION SUPERIORES A 1% DE HUMEDAD , EN ROJO LAS ZONAS DE SUBESTIMACION INFERIORES A 1% DE HUMEDAD Y EN COLOR CLARO LAS ZONAS CON PREDICCIONES ENTRE -1% Y 1% DE HUEMDAD

PROYECTAMOS LA PREDICCION A CORDENADAS GEOGRAFICAS

PREDICCION_P <- projectRaster(PREDICCION,crs="+init=epsg:4326")

PREDICCION_P

class       : RasterLayer 
dimensions  : 853, 912, 777936  (nrow, ncol, ncell)
resolution  : 0.00018, 0.000181  (x, y)
extent      : -75.85438, -75.69022, 4.417458, 4.571851  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : predicccion_SM 
values      : 76.31077, 79.72381  (min, max)

VISUALIZAMOS LA PREDICCION EN OPENSTREETMAP

pal <- colorNumeric(palette="YlGnBu", domain = c(75, 80),
  na.color = "transparent")
  
r <- PREDICCION_P

ra <- aggregate(r, fact=4, fun=mean)
leaflet() %>% addTiles() %>%
  addRasterImage(ra, colors = pal, opacity = 0.7) %>%
  addLegend(pal = pal, values = values(r),
    title = "PREDICCION DE HUMEDAD  (%) - MAYO 2015")

PREDICION EN OPENSTREETMAP