Extracción de datos SoilGrids

Ejemplo Nitrógeno (15-30 cm) Colombia

Agradecimientos

• Agradezco enormemente a Ivan Lizarazo por su documento guía publicado en Rpubs. Sus documentos son muy valiosos para todos aquellos que como yo queremos aprender un poco acerca de Sistemas de Información Geográfica con R.

SoilGrids

• SoilGrids

Bibliotecas

library(rgdal)
library(gdalUtils)
library(raster)
library(sf)
library(tidyverse)
library(ggspatial)

Raster Nitrógeno

• Datos SoilGrids
• ISRIC
• Artículo SoilGrids

nitro <- raster("https://files.isric.org/soilgrids/latest/data/nitrogen/nitrogen_15-30cm_mean.vrt")

nitro

class      : RasterLayer 
dimensions : 58034, 159246, 9241682364  (nrow, ncol, ncell)
resolution : 250, 250  (x, y)
extent     : -19949750, 19861750, -6147500, 8361000  (xmin, xmax, ymin, ymax)
crs        : +proj=igh +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source     : https://files.isric.org/soilgrids/latest/data/nitrogen/nitrogen_15-30cm_mean.vrt 
names      : nitrogen_15.30cm_mean 
values     : -32768, 32767  (min, max)

Shape Colombia

• Colombia Shape - Nivel geográfico departamentos

deptos <- st_read("data_shapes/MGN_DPTO_POLITICO.shp")

Reading layer `MGN_DPTO_POLITICO' from data source `D:\DocumentosEdimer\Github\Spatial-Data-Science\SoilGrids\data_shapes\MGN_DPTO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 33 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
Geodetic CRS:  WGS 84

• Proyección de coordenadas: transformación a proyección homolosena de Goode

igh <- '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
colombia_igh <- st_transform(deptos, igh)

colombia_igh

Simple feature collection with 33 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -9147200 ymin: -470815.4 xmax: -7442381 ymax: 1491094
CRS:           +proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs
First 10 features:
   DPTO_CCDGO DPTO_CNMBR DPTO_ANO_C                                     DPTO_ACT_A DPTO_NAREA
1          05  ANTIOQUIA       1886                  Constitucion Politica de 1886   62804.71
2          23    CÓRDOBA       1951              Ley 9 del 18 de Diciembre de 1951   25086.55
3          27      CHOCÓ       1947              Ley 13 del 3 de Noviembre de 1947   48353.22
4          70      SUCRE       1966                 Ley 47 del 8 de Agosto de 1966   10591.85
5          08  ATLÁNTICO       1910                                 Ley 21 de 1910    3313.81
6          13    BOLÍVAR       1886                  Constitucion Politica de 1886   26719.21
7          47  MAGDALENA       1964                                           1964   23135.94
8          20      CESAR       1967                    Ley 25  21 de junio de 1967   22565.31
9          44 LA GUAJIRA       1964 Acto Legislativo No. 1 de Diciembre 28 de 1964   20619.01
10         19      CAUCA       1857                            15 de junio de 1857   31242.91
   DPTO_CSMBL DPTO_VGNC Shape_Leng  Shape_Area                       geometry
1           3      2020  2397127.4 64104708675 MULTIPOLYGON (((-8537745 98...
2           3      2020  1084739.0 25771851296 MULTIPOLYGON (((-8487485 10...
3           3      2020  2305005.9 49095252907 MULTIPOLYGON (((-8641428 96...
4           3      2020   960645.8 10927904055 MULTIPOLYGON (((-8449360 11...
5           3      2020   286201.3  3452348633 MULTIPOLYGON (((-8384559 12...
6           3      2020  1820039.5 27531715587 MULTIPOLYGON (((-8426196 12...
7           3      2020  1213551.0 24044561103 MULTIPOLYGON (((-8298708 12...
8           3      2020  1410287.8 23352037403 MULTIPOLYGON (((-8229790 12...
9           3      2020  1211484.9 21584067313 MULTIPOLYGON (((-8050536 13...
10          3      2020  1553621.4 31435492565 MULTIPOLYGON (((-8515823 37...

• Área de interés (Colombia):

area <- st_bbox(colombia_igh)
area

      xmin       ymin       xmax       ymax 
-9147200.3  -470815.4 -7442381.4  1491094.3 

• Límites: formato numérico necesario para rgdal

ulx <- area$xmin
uly <- area$ymax
lrx <- area$xmax
lry <- area$ymin
limites <- c(ulx, uly, lrx, lry)
limites

      xmin       ymax       xmax       ymin 
-9147200.3  1491094.3 -7442381.4  -470815.4 

Nitrógeno Colombia (.tif)

• Obtención de datos con GDAL: este proceso se demora algunos minutos. Dará como resultado un archivo con formato .tif en la ruta especificada. Importante: es necesario tener instalado GDAL en nuestro computador (ver anexos para más información).

sg_url <- "/vsicurl/https://files.isric.org/soilgrids/latest/data/nitrogen/nitrogen_15-30cm_mean.vrt"
lfile <- "data_shapes/nitrogeno_colombia.tif"

gdal_translate(sg_url, lfile ,
               tr = c(250,250),
               projwin = limites,
               projwin_srs = igh,
               verbose = TRUE)

Checking gdal_installation...
Scanning for GDAL installations...
Checking Sys.which...
GDAL version 3.2.2
GDAL command being used: "C:\Program Files\GDAL\gdal_translate.exe" -tr 250 250 -projwin -9147200.27432083 1491094.27369704 -7442381.38468577 -470815.374159958 -of "GTiff" -projwin_srs "+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs" "/vsicurl/https://files.isric.org/soilgrids/latest/data/nitrogen/nitrogen_15-30cm_mean.vrt" "data_shapes/nitrogeno_colombia.tif"
Input file size is 159246, 580340...10...20...30ERROR 1: Request for 390313-406512 failed...40...50...60...70...80...90...100 - done.

NULL

Raster Nitrógeno Colombia

• Consultar unidades del nitrógeno

• Dividimos por 100 para tener las unidades originales (cg/kg)

nitro_colombia <- raster("data_shapes/nitrogeno_colombia.tif") / 100

plot(nitro_colombia)

Corte Colombia

• Corte con el shape de Colombia que tiene el mismo tipo de coordenadas:

corte_colombia <- nitro_colombia %>% 
  crop(colombia_igh) %>% 
  mask(colombia_igh)

plot(corte_colombia)

• Corte con shape obtenido con biblioteca raster:

# Mapa Colombia
ejemplo <- getData(name = "GADM", country = "COL", level = 0)

probando la URL 'https://biogeo.ucdavis.edu/data/gadm3.6/Rsp/gadm36_COL_0_sp.rds'
Content type 'à|ˆqþ' length 541789 bytes (529 KB)
downloaded 529 KB

ejemplo_sf <- st_as_sf(ejemplo)

# Reproyección
ejemplo_igh <- st_transform(ejemplo_sf, igh)

# Corte
corte_ejemplo <- nitro_colombia %>% 
  crop(ejemplo_igh) %>% 
  mask(ejemplo_igh)

plot(corte_ejemplo)

• Corte cambiando sistema de coordenadas:

nuevo_sistema <- "+proj=longlat +datum=WGS84 +no_defs"
prueba <- projectRaster(nitro_colombia, crs = nuevo_sistema)

# Corte
prueba_ejemplo <- prueba %>% 
  crop(ejemplo_sf) %>% 
  mask(ejemplo_sf)

plot(prueba_ejemplo)

Valle del Cauca

Filtro

valle <- colombia_igh %>% 
  filter(DPTO_CNMBR == "VALLE DEL CAUCA")

Corte

corte_valle <- nitro_colombia %>% 
  crop(valle) %>% 
  mask(valle)

plot(corte_valle, colNA = "black")

Mapa ggplot2

• Conversión a tibble: esta base de datos queda con 323589 observaciones.

valle_df <- corte_valle %>% 
  as("SpatialPixelsDataFrame") %>% 
  as_tibble()

valle_df %>% head()

ggplot() +
  geom_raster(data = valle_df, aes(x = x, y = y, fill = nitrogeno_colombia)) +
  geom_sf(data = valle, alpha = 0, fill = "transparent") +
  labs(fill = "Nitrógeno (cg/kg)",
       x = "Longitud", y = "Latitud") +
  annotation_north_arrow(location = "tl") +
  scale_fill_viridis_c(direction = -1) +
  theme_bw()

Antioquia

Filtro

antioquia <- colombia_igh %>% 
  filter(DPTO_CNMBR == "ANTIOQUIA")

Corte

corte_antioquia <- nitro_colombia %>% 
  crop(antioquia) %>% 
  mask(antioquia)

plot(corte_antioquia, colNA = "black")

Mapa ggplot2

• Conversión a tibble: esta base de datos queda con 995651 observaciones.

antioquia_df <- corte_antioquia %>% 
  as("SpatialPixelsDataFrame") %>% 
  as_tibble()

antioquia_df %>% head()

ggplot() +
  geom_raster(data = antioquia_df, aes(x = x, y = y, fill = nitrogeno_colombia)) +
  geom_sf(data = antioquia, alpha = 0) +
  labs(fill = "Nitrógeno (cg/kg)",
       x = "Longitud", y = "Latitud") +
  annotation_north_arrow(location = "tr") +
  scale_fill_gradientn(colours = rainbow(30)) +
  theme_bw()

Anexos

Recursos de información

• gdal.org
• gisinternals.com
• Tutorial instalación gdal Windows (inglés)
• Tutorial instalación gdal Windows (español)