Datos de humedad del suelo en Vichada

1. Introducción

Este cuaderno, ha sido creado desde RStudio Cloud e ilustra las funcionalidades de la biblioteca Smapr , en particular en este cuaderno se adquieren y se procesan los datos de humedad del suelo : SMAP (Soil Moisture Active-Passive) de la NASA. Además, se recortan los datos para que coincidan con el área de estudio y se generan visualizaciones simples de la humedad del suelo en el departamento del Vichada

2. Datos SMAP

La misión SMAP , mide la cantidad de agua en la superficie del suelo en todas partes de la Tierrae .Es decir, está diseñado para medir la humedad del suelo, cada 2-3 días

Algunos datos de esta misión SMAP , son proporcionados por el centro nacional de datos de nieve y hielo (NSIDC) . En la página del NSIDC se describen diferentes niveles de datos , principalmente se encuentran datos de nivel 3 que son productos compuestos globales y datos nivel 4 productos de datos globales modelados cada tres horas

Por otro lado,describen los productos de nivel 1 como mediciones de instrumentos brutos o calibrados y geolocalizados del radar y el radiometro SMAP, mientras que los datos de Nivel 2 , contienen recuperaciones del suelo derivadas de archivos auxiliares y los datos del nivel 1 .Estos 2 nieveles de datos tienen una resolución temporal de 49 minutos , que es el tiempo que tarda el satélite en completar media órbita de la tierra

Los productos del Nivel 3 son compuestos diarios de los datos de humedad del suelo del Nivel 2 y del estado de congelación/descongelación. Los productos del Nivel 4 proporcionan un modelo derivado de la humedad del suelo de la zona de las raíces y el intercambio de carbono neto del ecosistema.

Algunos de los datos que se obtienen del Satelite SMAP se listan a continuación

Productos SMAP

3.Autenticación en NSIDC

Para poder acceder los datos de la NASA a traves del NSDIC , se deben autenticar el nombre de usuario y la contraseña de Erathdata en el espacio de trabajo de RStudio

La forma más recomendada para la autenticación es usar Usar set_smap_credentials(‘nombre de usuario’, ‘contraseña’), una vez que se tenga instalada y cargada la libreria smapr.

4. Instalar y cargar paquetes

La isntalación y carga de paquetes se realizaron en la consola,los paquetes que fueron instalados fueron : Isoband (Desde el repositorio Git Hub),rhdf5 y smapr . Siendo esta última la libreria que permitirá obtener, descargar y extraer datos de nuestro interés

Autenticación

Como ya se habia mencionado , Se autenticará el acceso mediante el código que se mencionó en el numeral 3 , a traves de un usuario y contraseña de Earthdata . Se usa {r , include=FALSE} para que no se muestre este codigo en la salida

5. Encontrar datos SMAP

La función find_smap , devuelve los datos que estén disponibles según un criterio de búsqueda , que cosniste en la identificación del producto,la fecha y la versión a la que este pertenece , en este caso se buscarán los datos del Análisis global de la humedad en superficie y zona de raices

available_data <- find_smap(id = "SPL4SMAU", date = "2020-04-10", version = 4)
str(available_data)

'data.frame':   8 obs. of  3 variables:
 $ name: chr  "SMAP_L4_SM_aup_20200410T030000_Vv4030_001" "SMAP_L4_SM_aup_20200410T060000_Vv4030_001" "SMAP_L4_SM_aup_20200410T090000_Vv4030_001" "SMAP_L4_SM_aup_20200410T120000_Vv4030_001" ...
 $ date: Date, format: "2020-04-10" "2020-04-10" "2020-04-10" ...
 $ dir : chr  "SPL4SMAU.004/2020.04.10/" "SPL4SMAU.004/2020.04.10/" "SPL4SMAU.004/2020.04.10/" "SPL4SMAU.004/2020.04.10/" ...

6. Descargar Datos

En esta ocasión se descargara el producto mejorado de humedad del suelo de nivel 3 (L3) recuperadas por el radiómetro SMAP. Primero , se usará la función find_smap que se mencionó en el numeral anteriror

Archivos <- find_smap(id = "SPL3SMP_E", dates = "2020-04-28", version = 3)

Archivos

Luego, se descargan en la carpeta “suelo” que se crea en el directorio de trabajo

local_dir <- "./suelo"
downloads <- download_smap(Archivos[1, ], local_dir)

Downloading https://n5eil01u.ecs.nsidc.org/SMAP/SPL3SMP_E.003/2020.04.28/SMAP_L3_SM_P_E_20200428_R16515_001.h5
Downloading https://n5eil01u.ecs.nsidc.org/SMAP/SPL3SMP_E.003/2020.04.28/SMAP_L3_SM_P_E_20200428_R16515_001.qa
Downloading https://n5eil01u.ecs.nsidc.org/SMAP/SPL3SMP_E.003/2020.04.28/SMAP_L3_SM_P_E_20200428_R16515_001.h5.iso.xml

Se puede ver el contenido del objeto “downloads”

downloads

con la fución list_smap se puede ver el contenido (Formato HDF5) y los metadatos de lo que se descargó en el paso anterior

list_smap(downloads, all = TRUE)

$SMAP_L3_SM_P_E_20200428_R16515_001

NA

7. Visualiación de los Datos

Luego , podemos plotear la humedad del suelo de acuerdo al recorrido y los datos que recogió el satelite que se encuentra en el archivo descargado como “soil_mosture”

sm_raster1 <- extract_smap(downloads, "/Soil_Moisture_Retrieval_Data_AM/soil_moisture")
## sm_raster <- extract_smap("./soil/SMAP_L4_Antioquia/SMAP_L4_SM_lmc_00000000T000000_Vv4030_001", ## "/Geophysical_Data/sm_rootzone_wetness")
plot(sm_raster1, main = "Humedad del suelo nivel 3: 28 Marzo 2020")

Se puede realizar el mismo procedimiento pero para los datos existentes en otra hora del día

sm_raster <- extract_smap(downloads, "/Soil_Moisture_Retrieval_Data_PM/soil_moisture_pm")
## sm_raster <- extract_smap("./soil/SMAP_L4_Antioquia/SMAP_L4_SM_lmc_00000000T000000_Vv4030_001", ## "/Geophysical_Data/sm_rootzone_wetness")
plot(sm_raster, main = "Humedad del suelo nivel 3 : 28 Marzo 2020")

En los codigos anteriores , estos datos fueron almacenados en los objetos : sm_raster y sm_raster1 , se puede observar que hay en cada uno de esos objetos mediante los siguientes codigos:

sm_raster

class      : RasterLayer 
dimensions : 1624, 3856, 6262144  (nrow, ncol, ncell)
resolution : 9008.055, 9008.054  (x, y)
extent     : -17367530, 17367530, -7314540, 7314540  (xmin, xmax, ymin, ymax)
crs        : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /home/rstudio-user/.cache/smap/tmp.tif 
names      : SMAP_L3_SM_P_E_20200428_R16515_001 
values     : 0.01999992, 0.882247  (min, max)

sm_raster1

class      : RasterLayer 
dimensions : 1624, 3856, 6262144  (nrow, ncol, ncell)
resolution : 9008.055, 9008.054  (x, y)
extent     : -17367530, 17367530, -7314540, 7314540  (xmin, xmax, ymin, ymax)
crs        : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /home/rstudio-user/.cache/smap/tmp.tif 
names      : SMAP_L3_SM_P_E_20200428_R16515_001 
values     : 0.01999992, 0.8830189  (min, max)

Además , se puede indagar sobre el sistema de coordenadas de referencia de cada uno de los objetos

crs(sm_raster)

CRS arguments:
 +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0 

Posterior a esa operación, se lee el archivo shapefile correspondiente al departamento de vichada

Vichada <- read_sf("./ADMINISTRATIVO/MGN_DPTO_POLITICO.shp")
Vichada

Simple feature collection with 1 feature and 7 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -71.07793 ymin: 2.737109 xmax: -67.4098 ymax: 6.324317
CRS:            4326

¿Cuál es la extensión de los datos vectoriales?

extent(Vichada)

class      : Extent 
xmin       : -71.07793 
xmax       : -67.4098 
ymin       : 2.737109 
ymax       : 6.324317 

Los datos vectoriales están referidos a coordenadas geograficas , por tanto se deben reproyectar para que coincida con el sistema de referencia de los datos raster

cea_Vichada <- st_transform(Vichada, crs = crs(sm_raster))
cea_Vichada1 <- st_transform(Vichada, crs = crs(sm_raster1))

Se comprueba que ahora este objeto tiene el mismo sistema de referencia de los datos raster

crs(cea_Vichada)

[1] "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs "

crs(cea_Vichada1)

[1] "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs "

8.Recorte y Visualización de datos para el Área de estudio

Se hace uso de la función crop para que los datos raster de las diferentes horas del dia tengan la misma extensión del archivo shapefile

smap.sub <- crop(sm_raster, extent(cea_Vichada))
smap.sub1 <- crop(sm_raster1, extent(cea_Vichada1))

Se comprueba que los archivos tengan la extensión que se definió en el código anterior

smap.sub

class      : RasterLayer 
dimensions : 50, 39, 1950  (nrow, ncol, ncell)
resolution : 9008.055, 9008.054  (x, y)
extent     : -6855130, -6503816, 351314.1, 801716.8  (xmin, xmax, ymin, ymax)
crs        : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : SMAP_L3_SM_P_E_20200428_R16515_001 
values     : 0.08976348, 0.5132492  (min, max)

smap.sub1

class      : RasterLayer 
dimensions : 50, 39, 1950  (nrow, ncol, ncell)
resolution : 9008.055, 9008.054  (x, y)
extent     : -6855130, -6503816, 351314.1, 801716.8  (xmin, xmax, ymin, ymax)
crs        : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : SMAP_L3_SM_P_E_20200428_R16515_001 
values     : 0.08976348, 0.5132492  (min, max)

Se crea la función para convertir los valores en porcentajes,multiplicandolos x100

normalize <- function(x) {
  return(round(100*x))
}

Se aplica la función a los valores obtenidos

nuevo_smap = normalize(smap.sub)
nuevo_smap1 = normalize(smap.sub1)  

Se observa el contenido de los objetos

nuevo_smap

class      : RasterLayer 
dimensions : 50, 39, 1950  (nrow, ncol, ncell)
resolution : 9008.055, 9008.054  (x, y)
extent     : -6855130, -6503816, 351314.1, 801716.8  (xmin, xmax, ymin, ymax)
crs        : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 9, 51  (min, max)

nuevo_smap1

class      : RasterLayer 
dimensions : 50, 39, 1950  (nrow, ncol, ncell)
resolution : 9008.055, 9008.054  (x, y)
extent     : -6855130, -6503816, 351314.1, 801716.8  (xmin, xmax, ymin, ymax)
crs        : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 9, 51  (min, max)

Se crea un nuevo raster con la función “mask”

nuevo_smap2 <- mask(x = nuevo_smap, mask = cea_Vichada)
nuevo_smap1_2 <- mask(x = nuevo_smap1, mask = cea_Vichada1)

Se visualizan los datos con ayuda de la libreria leaflet, además se crea una nueva paleta de acuerdo a los valores que tome la humedad del suelo en el objeto “nuevo_smap2”

pal <- colorNumeric(c( "#FFFFCC", "#41B6C4","#0C2C84"), values(nuevo_smap2), na.color = "transparent")

leaflet() %>% addTiles() %>% addRasterImage(nuevo_smap2, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(nuevo_smap2), title = "Humedad del suelo en Vichada, 28 Marzo 2020 [%]")

Finalmente , se guardan los datos raster para ser usados en tareas posteriores

sm <- writeRaster(nuevo_smap2, filename="./suelo/sm_Vichada1_28_04_2020.tif", format="GTiff", datatype='INT1U', overwrite=TRUE)

Se guardan los datos para humedad del suelo

sm <- writeRaster(nuevo_smap1_2, filename="./suelo/sm_Vichada2_28_04_2020.tif", format="GTiff", datatype='INT1U', overwrite=TRUE)