Procesamiento de la base de datos global armonizada de suelos (HWSD) usando R
Ivan Lizarazo (Universidad Nacional de Colombia) y David Rossiter (Cornell University)
2017-05-29
Introducción
Esta nota explica cómo acceder y consultar la base de datos global armonizada de suelos (HWSD) [3] usando R, el programa estadístico de software libre [7]. Estas instrucciones permiten la integración de HWSD con otros datos geoespaciales y facilitan la ejecución de análisis estadísticos.
Esta nota describe cómo:
Acceder a los datos HWSD producidos por IIASA e importarlos a R;
Recuperar los datos HWSD para una ventana geográfica, usando una forma rectngular o poligonal;
Realizar la proyección cartográfica entre el sistema de referencia original, Plate Carrée, a un sistema de coordenadas plano como UTM;
Determinar el área cubierta por cada clase de suelo;
Guardar la ventana en el formato HWSD original y también como un raster proyectado;
Enlazar la base de datos de atributos con el raster y guardar los registros de la zona seleccionada en formato CSV o XLS;
Convertir la estructura raster original en estructura vectorial (polígonos);
Crear y desplegar loas atributos raster y los mapas vectoriales.
Existen cosas adicionales que se pueden hacer en R con los datos HWSD, incluyendo su integración con datos geoespaciales que están disponibles de manera gratuita, tales como modelos digitales de elevación e imágenes de satélite. Los lectores interesados pueden consultar el excelente texto de Bivand et al. [1] publicado en la serie de libros UseR! de Springer.
Todas las operaciones se realizan usando R, incluyendo el acceso directo a la base de datos HWSD cuyo formato original es MS-Access (extensión de archivo .mdb). Esto se realiza usando el paquete RODBC “R interface with Open Database Connectivity” [10].
Los procedimientos que se explican aquí usan diversos paquetes o librerías de R, incluyendo sp para el manejo de datos espaciales [1, 6], rgdal para la importación, la exportación y la transformación geométrica de datos espaciales [4], raster para trabajar con imágenes y datos en formato raster (grid) [2], y, como ya se indicó, RODBC para el manejo de la base de datos [10].
Estos paquetes deben ser cargados antes de usarlos por primera vez, tal como se muestra en el código.
Nota: El código incluido en este documento fue probado con la version 3.3.3 de R (2017-03-06) y con paquetes de dicha versión corriendo en Windows 10 Pro. El texto y la salida gráficas fueron escritas como un archivo Rmd usando el programa RStudio 1.0.136 que fue compilado mediante Knit para obtener este archivo PDF.
Si usted tiene sugerencias u observaciones sobre esta nota puede escribir a Ivan Lizarazo y/o a David G. Rossiter.
Esta nota es una traducción al español de la nota técnica “Processing the Harmonized World Soil Database (Version 1.2) in R” escrita por D. G. Rossiter (2014) ©.
La nota técnica original ha sido adaptada para trabajar directamente con la base de datos en formato MS-Access, en lugar del formato SQLite usado en dicha nota. Como los datos HWSD son proporcionados por IIASA en formato .mdb, este cambio puede ser útil. Adicionalmente, en esta adaptación los ejemplos se realizan en una zona geográfica localizada al oriente de Bogotá, en Colombia, América del Sur, en lugar de la zona geográfica utilizada en la nota original, localizada en las provincias de Jiangsu y Anhui en China.
1. Importación de datos HWSD en R
Tarea 1 : Descargar la base de datos HWSD de la página web de IIASA.
Los datos HWSD se pueden descargar desde IIASA. Esos datos comprenden tres archivos:
HWSD_raster.zip Datos en formato raster (.bil, .blw, .hdr)
HWSD.mdb Soil Base de datos de atributos del suelo en formato MS Access (.mdb)
HWSD_META.mdb Metadatos de atributos del suelo en formato MS Access (.mdb)
Tarea 2 : Descomprimir el archivo HWSD_RASTER.zip.
Al finalizar esta tarea se obtiene un subdirectorio HWSD_RASTER con tres archivos: un archivo raster con datos en formato intercalado por banda (“band interleaved”) (hwsd.bil, 1.7 Gb), un archivo pequeño con la extensión y la resolución de los datos (hwsd.blw), y un archivo con los datos de encabezado (hwsd.hdr). Estos últimos dos archivos son consultados automáticamente cuando se realiza la importación de los datos.
Tarea 3 : Importar los datos raster a R.
El paquete raster permite trabajar con archivos raster de gran tamaño, como el que nos ocupa, debido a que solamente guarda en memoria los datos que necesita, mientras que el resto permanece en el disco. La función raster asocia un objeto de R con el archivo que se encuentra en disco. El formato .bil es uno de los varios formatos que se pueden importar usando el comando raster. El paquete raster depende del paquete sp que se carga automáticamente si es requerido. El paquete raster se carga usando la función require que habilita el uso de cualquier paquete que no esté en el espacio de trabajo. Antes de habilitar los paquetes, se debe definir el directorio de trabajo. Para ello se usará el comando setwd(ruta_a_su_directorio). Tenga en cuenta que, en Windows, la ruta de un directorio se indica con // en lugar de /.
setwd("C://Users//Iván Lizarazo//Documents//UN//datos//hwsd")
require(sp)## Loading required package: sprequire(raster)## Loading required package: rasterhwsd <- raster("./hwsd.bil")
hwsd## class : RasterLayer
## dimensions : 21600, 43200, 933120000 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : C:\Users\Iván Lizarazo\Documents\UN\datos\hwsd\hwsd.bil
## names : hwsd
## values : 0, 65535 (min, max)Tarea 4: Examinar las propiedades del archivo raster.
El paquete raster ofrece comandos útiles para esta tarea, cuyos nombres son obvios:
ncol(hwsd)## [1] 43200nrow(hwsd)## [1] 21600res(hwsd)## [1] 0.008333333 0.008333333extent(hwsd)## class : Extent
## xmin : -180
## xmax : 180
## ymin : -90
## ymax : 90projection(hwsd)## [1] NAComo se observa, el archivo raster no tiene información de su referencia espacial, por eso la respuesta del comando projection es NA (“not available”). Si examinamos la documentación de HWSD [3] podemos conocer que los datos usan la proyección Plate Carrée y el datum WGS84; es decir que los datos de latitud y longitud se trasladan de manera directa a una grilla, de manera que la geometría obtenida tiene mayor distorsión a medida que las posiciones están más cerca de los polos.
Tarea 5: Asociar la información de referencia espacial a los datos raster.
Para expresar la información de proyección se puede usar la función proj4string que utiliza la sintaxis del formato PROJ4 [5]. Podemos indicar que la “proyección” no existe, es decir que los datos están en coordenadas geográficas. Adicionalmente, indicamos el datum, el elipsoide y los parámetros de transformación a WGS84. Note que el nombre el datum WGS84 debe estar en mayúsculas:
require(rgdal)## Loading required package: rgdal## rgdal: version: 1.2-6, (SVN revision 651)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: c:/Rlibs/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: c:/Rlibs/rgdal/proj
## Linking to sp version: 1.2-4(proj4string(hwsd) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")## [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"2.Selección de la región de interés
Como la base de datos es muy grande, usualmente trabajamos en una región geográfica específica.
2.1 Selección mediante rectángulo envolvente
El paquete raster permite recortar un archivo raster usando una zona de interés (“extent”). Esta zona puede ser extraída a partir del rectángulo envolvente de cualquier objeto sp. También puede ser especificado de manera directa usando la función extent. En este nota vamos a usar esa función para seleccionar una zona de 3° by 3° con centro cerca de Bogotá, Colombia. Luego, vamos a recortar el archivo raster usando la función crop:
hwsd.crop <- crop(hwsd, extent(c(-76, -73, 4, 7)))
nrow(hwsd.crop)## [1] 360ncol(hwsd.crop)## [1] 360bbox(hwsd.crop)## min max
## s1 -76 -73
## s2 4 7La función unique muestra los códigos que identifican las clases de suelo existentes en la zona:
unique(hwsd.crop)## [1] 16002 16003 16009 16010 16013 16014 16015 16016 16017 16018 16024
## [12] 16028 16030 16031 16039 16042 16047 16049 16053 16054Este es el único contenido de un archivo raster: cada celda tiene un código que permite enlazarse con la base de datos, tal como se indica más adelante.
Tarea 6: Despliegue de los datos raster con una paleta de color apropiada.
Puede observarse que existen 20 clases de suelos. Este número es relativamente alto para representar con diferentes colores. Una manera de hacerlo es utilizar una paleta con una rampa colores continuos:
plot(hwsd.crop, main = "Unidades de suelos HWSD", col = bpy.colors(length(unique(hwsd.crop))))
La representación obtenida parece buena, toda vez que los códigos parecen estar ordenados por similaridad de suelos. También podemos los primeros dos dígitos de los códigos de suelos, lo cual es presumiblemente apropiado. Para remover las “centenas” podemos usar el operador de división entre enteros %/%. El paquete RColorBrewer proporciona diversas paletas de color; aquí se seleccionó la paleta de nombre “Accent” que realza las diferencias entre clases. Para ello usamos la función brewer.pal:
hwsd.crop3 <- (hwsd.crop%/%10)
freq(hwsd.crop3)## value count
## [1,] 1600 13972
## [2,] 1601 56473
## [3,] 1602 3267
## [4,] 1603 12374
## [5,] 1604 40510
## [6,] 1605 3004require(RColorBrewer)## Loading required package: RColorBrewerplot(hwsd.crop3, main = "Unidades de suelos HWSD", col = brewer.pal(length(unique(hwsd.crop3)),
"Accent"))
Otra opción de representación es el uso del paquete leaflet:
library(raster)
library(leaflet)
pal <- colorNumeric(palette = "YlGnBu", domain = unique(hwsd.crop),
na.color = "transparent")
leaflet() %>% addTiles() %>% addRasterImage(hwsd.crop, colors = pal,
opacity = 0.5) %>% addLegend(pal = pal, values = unique(hwsd.crop),
title = "Unidades de suelos HWSD")En cualquier caso, los datos en coordenadas geográficas pueden aparecer distorsionados. Podemos observar el efecto de la proyección usando el método projectRaster y especificando un sistema de referencia de coordenadas (CRC) de salida. Observe que vamos a usar el método de remuestro del vecino más cercano (method=“ngb”) debido a que los datos corresponden a una clasificación de suelos.
Nuestro sistema de coordenadas planas será MAGNA-SIRGAS Colombia Bogota Zone, teniendo en cuenta que ese es el sistema oficial de referencia en Colombia. En la página Spatial Reference podemos ver que el código EPSG de dicho sistema es 3116. Allí mismo podemos obtener la cadena de texto proj4 correspondiente a este sistema:
print(paste("UTM zone:", utm.zone <- floor((sum(bbox(hwsd.crop)[1,
])/2 + 180)/6) + 1))## [1] "UTM zone: 18"proj4string.3116 <- "+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"
hwsd.plane <- projectRaster(hwsd.crop, crs = proj4string.3116,
method = "ngb")
hwsd.plane## class : RasterLayer
## dimensions : 370, 371, 137270 (nrow, ncol, ncell)
## resolution : 923, 922 (x, y)
## extent : 781852.1, 1124285, 929730.6, 1270871 (xmin, xmax, ymin, ymax)
## coord. ref. : +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
## data source : in memory
## names : hwsd
## values : 16002, 16054 (min, max)plot(hwsd.plane, main = "Unidades de suelos HWSD", col = brewer.pal(11,
"Accent"), asp = 1)## Warning in brewer.pal(11, "Accent"): n too large, allowed maximum for palette Accent is 8
## Returning the palette you asked for with that many colorsgrid()
Ahora que los datos están en una representación plana, geométricamente correcta, podemos calcular el área cubierta por cada clase de suelo, lo mismo que el área total de la zona de estudio, en km2:
(cell.dim <- res(hwsd.plane))## [1] 923 922(cell.area <- cell.dim[1] * cell.dim[2]/10^4)## [1] 85.1006(tmp <- cbind(freq(hwsd.plane)[, 1], freq(hwsd.plane)[, 2] *
cell.area/10^2))## [,1] [,2]
## [1,] 16002 8736.4276
## [2,] 16003 2309.6303
## [3,] 16009 844.1980
## [4,] 16010 4289.0702
## [5,] 16013 6783.3688
## [6,] 16014 10435.8866
## [7,] 16015 13070.6012
## [8,] 16016 559.1109
## [9,] 16017 6586.7864
## [10,] 16018 6363.8229
## [11,] 16024 1338.6324
## [12,] 16028 1432.2431
## [13,] 16030 2263.6760
## [14,] 16031 2131.7700
## [15,] 16039 6158.7304
## [16,] 16042 167.6482
## [17,] 16047 28708.6874
## [18,] 16049 5646.4248
## [19,] 16053 188.0723
## [20,] 16054 2368.3497
## [21,] NA 6434.4564ix <- which(is.na(tmp[, 1]))
sum(tmp[-ix, 2])## [1] 110383.1rm(cell.dim, cell.area, tmp, ix)Observe que el área de cada una de las celdas equivale aproximadamente a 85 ha. La zona de interés cubre alrededor de 110 000 km2. Observe también que existen algunas celdas con valor NA. Estas celdas corresponden a los bordes del raster proyectado, las cuales se requieren para asegurar que la zona mantenga su geometría rectangular.
Como ya revisamos el raster proyectado, podemos removerlo:
rm(hwsd.plane)Volviendo al raster no proyectado, podemos consultar la clase en cualquier sitio usando la función click. Luego de llamar esta función, haga clic con el mouse en una celda en la imagen desplegada. Ello regresará las coordenadas del punto, y, si el argumento opcional id es definido como verdadero, el código de la celda será devuelto. Por ejemplo, al hacer clic en el oriente de Bogotá, dentro del Parque Nacional Chingaza (Latitud = 4.6875, Longitud =-73.87083) se obtendrán los siguientes valores:
plot(hwsd.crop, main = "Unidades de Suelos HWSD", col = bpy.colors(length(unique(hwsd.crop))))
# xy <- click(hwsd.crop, n = 1, id = TRUE, xy = TRUE, type =
# 'p') print(xy)
xy <- data.frame(x = -73.87083, y = 4.6875, value = 16017)
hwsd.id <- xy["value"]
hwsd.id## value
## 1 16017Observe que el código que llama la función click está comentariado. Esto se hace así para facilitar la compilación de este documento. Note que las coordenadas están en grados decimales y que el resultado es un identificador o código que indica la unidad cartográfica de suelos correspondiente a la celda de interés.
2.2 Selección mediante polígono envolvente
Otra manera de seleccionar un subconjunto de la base de datos utiliza los límites poligonales de una región, por ejemplo, un departamento, que es una entidad política-administrativa que agrupa diversos municipios.
Tarea 7: Crear un objeto SpatialPolygons a partir de los límites del Departamento del Meta
Podemos obtener un archivo con los límites del departamento del Meta en la página SIGOT del Instituto Geográfico Agustín Codazzi (IGAC). Para leer el shapefile correspondiente a este departamento usamos la función shapefile del paquete raster. El objeto resultante es un objeto de la clase SpatialPolygonsDataFrame. Este objeto se puede convertir a un objeto de la clase SpatialPolygons.
meta <- shapefile("meta.shp")
meta## class : SpatialPolygonsDataFrame
## features : 1
## extent : 1110999, 1470026, 965835.4, 1190167 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200417 +lon_0=-74.07750791700001 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs
## variables : 1
## names : NOM_DEPART
## min values : CASANARE
## max values : CASANAREmeta.boundary <- as(meta, "SpatialPolygons")
meta.boundary## class : SpatialPolygons
## features : 1
## extent : 1110999, 1470026, 965835.4, 1190167 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200417 +lon_0=-74.07750791700001 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defsbbox(meta.boundary)## min max
## x 1110999.4 1470026
## y 965835.4 1190167Tarea 8 : Extraer la porción de la base de datos HWSD que corresponde al polígono de interés
Cuando se trabaja con objetos RasterLayer solamente podemos extraer áreas rectangulares. Por ello, lo primero que vamos a hacer es extraer el rectángulo envolvente de Meta para extraer ese departamento y algunas áreas de los departamentos vecinos a un objeto SpatialPixelDataFrame, tal como lo requiere el paquete sp. Enseguida, vamos a usar la función over para extraer solamente los pixeles que están dentro del departamento del Meta. Esta función retorna NA para las celdas que están por fuera del departamento, y el identificador ID del polígono para las celdas que están dentro. Enseguida, seleccionamos las celdas que no son NA, es decir, aquellas que tienen un identificador.
# crop Raster* with Spatial* object
proj_hwsd <- proj4string(hwsd)
crs_hwsd <- CRS(proj_hwsd)
meta.bound <- spTransform(meta.boundary, crs_hwsd)
hwsd.meta.box <- crop(hwsd, bbox(meta.bound))
hwsd.meta.box.sp <- as(hwsd.meta.box, "SpatialPixelsDataFrame")
sort(unique(hwsd.meta.box.sp$hwsd))## [1] 16014 16015 16017 16018 16024 16025 16026 16028 16029 16034 16039
## [12] 16047 16049 16051 16054ix <- over(hwsd.meta.box.sp, meta.bound)
hwsd.meta.sp <- hwsd.meta.box.sp[!is.na(ix), ]
(meta.id <- sort(unique(hwsd.meta.sp$hwsd)))## [1] 16017 16018 16025 16026 16028 16029 16047 16051 16054spplot(hwsd.meta.box.sp, main = "Suelos HWSD, Meta - Rectángulo envolvente",
col.regions = topo.colors(64), scales = list(draw = TRUE))
spplot(hwsd.meta.sp, main = "Suelos HWSD, Meta - Polígono envolvente",
col.regions = topo.colors(64), scales = list(draw = TRUE))
Los datos pueden, entonces, ser convertidos nuevamente en un objeto RasterLayer:
hwsd.meta <- as(hwsd.meta.sp, "RasterLayer")
rm(hwsd.meta.box, hwsd.meta.box.sp, ix)3. Base de datos de atributos
Enseguida vamos a usar el paquete RODBC [10] para leer la base de datos MS Access databases (extensión .mdb). Tenga en cuenta que este paquete solamente funciona con Windows 32 bit. Por consiguiente, si su sistema es Windows 64 bits, usted debe presionar la tecla Control durante el inicio de RStudio con el objeto de desplegar la ventana de seleccion de la versión R. De esta manera, podrá seleccionar la versión de R para 32 bits.
Las instrucciones para usar el paquete RODBC que se usan en esta nota, se basan en [10] así como en el código descrito en R-bloggers.
Tarea 9: Conectarse a la base de datos de atributos HWSD and listar las tablas.
Primero hay que cargar el paquete RODBC. Enseguida usamos la función odbcConnectAccess para establecer la conexión con la base de datos de interés. La función sqlTables lista las tablas relacionales que existen en la base de datos.
# Load RODBC package
library(RODBC)
# Connect to Access db
tcon <- odbcConnectAccess("C://Users//Iván Lizarazo/Documents//UN//datos//hwsd/HWSD.mdb")
# List tables
tbls <- sqlTables(tcon)
tbls$TABLE_NAME## [1] "MSysAccessObjects" "MSysAccessXML"
## [3] "MSysACEs" "MSysIMEXColumns"
## [5] "MSysIMEXSpecs" "MSysNavPaneGroupCategories"
## [7] "MSysNavPaneGroups" "MSysNavPaneGroupToObjects"
## [9] "MSysNavPaneObjectIDs" "MSysObjects"
## [11] "MSysQueries" "MSysRelationships"
## [13] "D_ADD_PROP" "D_AWC"
## [15] "D_COVERAGE" "D_DRAINAGE"
## [17] "D_IL" "D_ISSOIL"
## [19] "D_PHASE" "D_ROOTS"
## [21] "D_SWR" "D_SYMBOL"
## [23] "D_SYMBOL74" "D_SYMBOL85"
## [25] "D_SYMBOL90" "D_TEXTURE"
## [27] "D_USDA_TEX_CLASS" "HWSD_DATA"
## [29] "HWSD_SMU" "NTEMPO10"
## [31] "OVENTANA" "TEMPO10"
## [33] "VENTANA" "HWSD_Q"# sqlDrop(tcon, 'VENTANA') sqlDrop(tcon, 'NTEMPO')Observe que la base de datos HWSD tiene 18 archivos. La tabla HWSD_DATA contiene la lista de las unidades de suelo, mientras que el resto de tablas son tablas de referencia (“lookup tables”) de los códigos de los atributos.
Tarea 10: Desplegar la estructura de la tabla principal.
La función sqlColumns() regresa los nombres de las columnas de una tabla, es decir, los nombres de las variables. Esa función requiere dos argumentos: el nombre de la conexión y el nombre de la tabla. Los siguientes comandos regresan el nombre de las columnas y el tipo de datos correspondientes a la tabla HWSD_DATA.
sqlColumns(tcon, "HWSD_DATA")$COLUMN_NAME## [1] "ID" "MU_GLOBAL" "MU_SOURCE1"
## [4] "MU_SOURCE2" "ISSOIL" "SHARE"
## [7] "SEQ" "SU_SYM74" "SU_CODE74"
## [10] "SU_SYM85" "SU_CODE85" "SU_SYM90"
## [13] "SU_CODE90" "T_TEXTURE" "DRAINAGE"
## [16] "REF_DEPTH" "AWC_CLASS" "PHASE1"
## [19] "PHASE2" "ROOTS" "IL"
## [22] "SWR" "ADD_PROP" "T_GRAVEL"
## [25] "T_SAND" "T_SILT" "T_CLAY"
## [28] "T_USDA_TEX_CLASS" "T_REF_BULK_DENSITY" "T_BULK_DENSITY"
## [31] "T_OC" "T_PH_H2O" "T_CEC_CLAY"
## [34] "T_CEC_SOIL" "T_BS" "T_TEB"
## [37] "T_CACO3" "T_CASO4" "T_ESP"
## [40] "T_ECE" "S_GRAVEL" "S_SAND"
## [43] "S_SILT" "S_CLAY" "S_USDA_TEX_CLASS"
## [46] "S_REF_BULK_DENSITY" "S_BULK_DENSITY" "S_OC"
## [49] "S_PH_H2O" "S_CEC_CLAY" "S_CEC_SOIL"
## [52] "S_BS" "S_TEB" "S_CACO3"
## [55] "S_CASO4" "S_ESP" "S_ECE"sqlColumns(tcon, "HWSD_DATA")$TYPE_NAME## [1] "INTEGER" "INTEGER" "VARCHAR" "INTEGER" "BYTE" "REAL"
## [7] "BYTE" "VARCHAR" "SMALLINT" "VARCHAR" "SMALLINT" "VARCHAR"
## [13] "SMALLINT" "BYTE" "SMALLINT" "SMALLINT" "BYTE" "BYTE"
## [19] "BYTE" "BYTE" "BYTE" "BYTE" "BYTE" "SMALLINT"
## [25] "SMALLINT" "SMALLINT" "SMALLINT" "BYTE" "DOUBLE" "DOUBLE"
## [31] "DOUBLE" "DOUBLE" "DOUBLE" "DOUBLE" "DOUBLE" "DOUBLE"
## [37] "DOUBLE" "DOUBLE" "DOUBLE" "DOUBLE" "SMALLINT" "SMALLINT"
## [43] "SMALLINT" "SMALLINT" "BYTE" "DOUBLE" "DOUBLE" "DOUBLE"
## [49] "DOUBLE" "DOUBLE" "DOUBLE" "DOUBLE" "DOUBLE" "DOUBLE"
## [55] "DOUBLE" "DOUBLE" "DOUBLE"Tarea 11 : Determinar el número de registros en la tabla principal
Podemos usar la función SQL count para contar el número de registros seleccionados (en este caso, todos los registros se representan con un asterisco) y asignar un nombre al usando el operador as. La selección de todos los registros se puede realizar omitiendo el uso de la cláusula where.
sqlQuery(tcon, "select count(*) as grid_total from HWSD_DATA")## grid_total
## 1 48148Tarea 12: Desplegar el ID, el código de la unidad de suelos, indicar si es una unidad de suelo o no, el porcentaje de la unidad, el código de la clase FAO 1990 y los código de textura de la capa superior del suelo, para los primeros diez registros de la tabla principal.
(display.fields <- c("ID", "MU_GLOBAL", "ISSOIL", "SHARE", "SU_CODE90",
"SU_SYM90", "T_USDA_TEX_CLASS"))## [1] "ID" "MU_GLOBAL" "ISSOIL"
## [4] "SHARE" "SU_CODE90" "SU_SYM90"
## [7] "T_USDA_TEX_CLASS"tmp <- sqlQuery(tcon, paste("select", paste(display.fields, collapse = ", "),
"from HWSD_DATA"))
dim(tmp)## [1] 48148 7print(tmp[1:10, display.fields])## ID MU_GLOBAL ISSOIL SHARE SU_CODE90 SU_SYM90 T_USDA_TEX_CLASS
## 1 2 7002 0 100 202 HD NA
## 2 3 7003 0 100 198 WR NA
## 3 4 7004 0 100 89 HSf 3
## 4 5 7005 0 100 199 GG NA
## 5 6 7006 1 70 35 ANz 11
## 6 7 7006 1 20 32 ANh 11
## 7 8 7006 1 10 37 ANi 9
## 8 9 7007 1 80 35 ANz 11
## 9 10 7007 1 20 32 ANh 11
## 10 11 7008 1 90 35 ANz 11rm(tmp)Podemos observar que algunas unidades (p.ej. 7001) no corresponden a suelos. Otras unidades (p.ej. 7004) tienen solamente un componente y, otros (p.ej. 7006) tienen varios componentes, con sus respectivas proporciones.
Tarea 13: Desplegar la estructura de la tabla de referencia de las clases de suelo FAO 1990.
En la documentación de los datos HWSD se indica que las tablas de referencia tienen nombres que usan la convención D_*. La tabla que describe las clases FAO 1990 se llama D_SYMBOL90. Como esta table es muy pequeña podemos leerla en memoria seleccionando todas las filas. Luego, podemos examinar su estructura.
str(sqlQuery(tcon, "select * from D_SYMBOL90"))## 'data.frame': 193 obs. of 3 variables:
## $ CODE : int 1 2 3 4 5 6 7 8 9 10 ...
## $ VALUE : Factor w/ 193 levels "ACRISOLS","Albic Arenosols",..: 64 45 14 36 141 181 177 165 92 46 ...
## $ SYMBOL: Factor w/ 193 levels "AC","ACf","ACg",..: 55 58 56 57 59 62 61 60 72 75 ...Tarea 14: Mostrar el registro correspondiente a la celda identificada en la sección anterior.
Nuevamente, podemos usar la función sqlQuery, pero ahora usando una cadena de texto para construir la consulta que nos permitirá encontrar el registro deseado. Observe que usamos la función paste para crear la cadena de texto usando palabras fijas (entre comillas) lo mismo que texto obtenido a partir de una variable. En este caso, el código de la unidad de suelos que se guardó en la variable hwsd.id durante la consulta interactiva del mapa que se realizó anteriormente.
(tuple <- sqlQuery(tcon, paste("select * from HWSD_DATA where MU_GLOBAL = ",
hwsd.id)))## ID MU_GLOBAL MU_SOURCE1 MU_SOURCE2 ISSOIL SHARE SEQ SU_SYM74
## 1 15632 16017 CO26 NA 1 60 1 NA
## 2 15633 16017 CO26 NA 1 40 2 NA
## SU_CODE74 SU_SYM85 SU_CODE85 SU_SYM90 SU_CODE90 T_TEXTURE DRAINAGE
## 1 NA NA NA LPu 102 2 3
## 2 NA NA NA ANu 34 2 4
## REF_DEPTH AWC_CLASS PHASE1 PHASE2 ROOTS IL SWR ADD_PROP T_GRAVEL T_SAND
## 1 30 5 NA NA NA NA NA 0 1 60
## 2 100 1 NA NA NA NA NA 0 0 68
## T_SILT T_CLAY T_USDA_TEX_CLASS T_REF_BULK_DENSITY T_BULK_DENSITY T_OC
## 1 29 11 11 1.53 1.10 6.62
## 2 26 6 11 1.63 0.52 5.49
## T_PH_H2O T_CEC_CLAY T_CEC_SOIL T_BS T_TEB T_CACO3 T_CASO4 T_ESP T_ECE
## 1 4.1 40 28 38 10.6 0 0 4 0
## 2 5.0 140 29 14 4.1 0 0 1 0
## S_GRAVEL S_SAND S_SILT S_CLAY S_USDA_TEX_CLASS S_REF_BULK_DENSITY
## 1 NA NA NA NA NA NA
## 2 0 68 28 4 11 1.69
## S_BULK_DENSITY S_OC S_PH_H2O S_CEC_CLAY S_CEC_SOIL S_BS S_TEB S_CACO3
## 1 NA NA NA NA NA NA NA NA
## 2 0.62 4.1 5.4 140 29 19 5.5 0
## S_CASO4 S_ESP S_ECE
## 1 NA NA NA
## 2 0 1 0tuple$SU_SYM90## [1] LPu ANu
## Levels: ANu LPuEl resultado anterior muestra el código. Podemos encontrar el nombre correspondiente en la tabla de referencia:
sqlQuery(tcon, paste("select * from D_SYMBOL90 where symbol='",
tuple$SU_SYM90, "'", sep = ""))## CODE VALUE SYMBOL
## 1 102 Umbric Leptosols LPuUn agrólogo que conozca el área de estudio, nos puede confirmar si la clase obtenida es razonable o no.
Ahora podemos elaborar un mapa de propiedades en la ventana de interés.
Tarea 15: Extraer las unidades de suelos en la región.
Una manera de extraer los registros requeridos de la base de datos es crear una tabla con la lista de las unidades cartográficas existentes en la zona, y luego, usar esa tabla como criterio de selección usando un JOIN. La función sqlSave crea una tabla, para ello se requiere un data frame de R como valor inicial. A partir de ese data frame, se infiere la estructura de la tabla.
## sqlDrop(tcon, 'NVENTANA')
sqlSave(channel = tcon, dat = data.frame(SMU_ID = unique(hwsd.crop)),
tablename = "NVENTANA", addPK = TRUE)
qry <- "SELECT * FROM NVENTANA"
sqlColumns(tcon, "NVENTANA")$COLUMN_NAME## [1] "rownames" "SMU_ID"datos <- sqlQuery(tcon, qry)
datos## rownames SMU_ID
## 1 1 16002
## 2 2 16003
## 3 3 16009
## 4 4 16010
## 5 5 16013
## 6 6 16014
## 7 7 16015
## 8 8 16016
## 9 9 16017
## 10 10 16018
## 11 11 16024
## 12 12 16028
## 13 13 16030
## 14 14 16031
## 15 15 16039
## 16 16 16042
## 17 17 16047
## 18 18 16049
## 19 19 16053
## 20 20 16054Enseguida, vamos a realizar la unión usando el atributo común. La nueva tabla no aporta ningún atributo nuevo. También vamos a mostrar cómo ordenar los resultados usando el símbolo de unidades de FAO 1990.
Nota: La cláusula T.* selecciona los atributos de la tabla HWSD_DATA; ello se representa por T en la instrucción de unión. No se requiere los atributos de la tabla con la lista de las unidades de suelo que existen en la región, toda vez que la tabla HWSD_DATA tiene dichos códigos en el atributo MU_GLOBAL.
records <- sqlQuery(tcon, "SELECT * FROM HWSD_DATA AS T
INNER JOIN NVENTANA
AS U ON T.MU_GLOBAL = U.SMU_ID")
dim(records)## [1] 41 59head(records)## ID MU_GLOBAL MU_SOURCE1 MU_SOURCE2 ISSOIL SHARE SEQ SU_SYM74
## 1 15606 16002 CO12 NA 1 55 1 NA
## 2 15607 16002 CO12 NA 1 45 2 NA
## 3 15608 16003 CO13 NA 1 60 1 NA
## 4 15609 16003 CO13 NA 1 40 2 NA
## 5 15618 16009 CO19 NA 1 100 1 NA
## 6 15619 16010 CO2 NA 1 40 1 NA
## SU_CODE74 SU_SYM85 SU_CODE85 SU_SYM90 SU_CODE90 T_TEXTURE DRAINAGE
## 1 NA NA NA CMd 63 2 4
## 2 NA NA NA CMo 68 3 4
## 3 NA NA NA CMe 62 1 5
## 4 NA NA NA LPe 98 1 4
## 5 NA NA NA ANu 34 2 4
## 6 NA NA NA GLe 10 3 1
## REF_DEPTH AWC_CLASS PHASE1 PHASE2 ROOTS IL SWR ADD_PROP T_GRAVEL T_SAND
## 1 100 1 NA NA NA NA NA 0 0 29
## 2 100 1 NA NA NA NA NA 0 3 39
## 3 100 3 NA NA NA NA NA 0 1 75
## 4 30 5 NA NA NA NA NA 0 19 87
## 5 100 1 NA NA NA NA NA 0 0 45
## 6 100 1 NA NA NA NA NA 0 0 13
## T_SILT T_CLAY T_USDA_TEX_CLASS T_REF_BULK_DENSITY T_BULK_DENSITY T_OC
## 1 37 34 5 1.30 1.31 0.45
## 2 23 38 5 1.31 1.27 0.37
## 3 13 12 11 1.54 1.50 0.61
## 4 8 5 12 1.70 1.70 1.01
## 5 34 21 9 1.40 0.65 4.36
## 6 33 54 3 1.20 1.26 0.99
## T_PH_H2O T_CEC_CLAY T_CEC_SOIL T_BS T_TEB T_CACO3 T_CASO4 T_ESP T_ECE
## 1 4.5 22 9 28 2.5 0.0 0 1 0.0
## 2 4.6 19 8 38 3.0 0.0 0 2 0.0
## 3 6.2 45 7 87 6.1 0.0 0 2 0.0
## 4 6.3 78 7 80 5.6 0.1 0 3 0.1
## 5 5.8 33 21 14 2.9 0.0 0 1 0.0
## 6 5.5 58 36 65 23.4 0.0 0 2 0.0
## S_GRAVEL S_SAND S_SILT S_CLAY S_USDA_TEX_CLASS S_REF_BULK_DENSITY
## 1 0 32 34 34 5 1.31
## 2 2 25 27 48 3 1.24
## 3 8 78 12 10 11 1.58
## 4 NA NA NA NA NA NA
## 5 0 31 37 32 5 1.32
## 6 0 7 47 46 2 1.22
## S_BULK_DENSITY S_OC S_PH_H2O S_CEC_CLAY S_CEC_SOIL S_BS S_TEB S_CACO3
## 1 1.30 0.40 4.6 21 9 29 2.6 0.0
## 2 1.14 0.10 4.8 23 11 45 5.0 0.0
## 3 1.63 0.21 6.5 65 7 87 6.1 0.2
## 4 NA NA NA NA NA NA NA NA
## 5 0.81 1.88 5.9 50 20 38 7.6 0.0
## 6 1.48 0.20 6.7 53 26 88 22.9 0.3
## S_CASO4 S_ESP S_ECE rownames SMU_ID
## 1 0.0 1 0.0 1 16002
## 2 0.0 2 0.0 1 16002
## 3 0.1 3 0.1 2 16003
## 4 NA NA NA 2 16003
## 5 0.0 1 0.0 3 16009
## 6 0.1 2 0.1 4 16010head(records)[, display.fields]## ID MU_GLOBAL ISSOIL SHARE SU_CODE90 SU_SYM90 T_USDA_TEX_CLASS
## 1 15606 16002 1 55 63 CMd 5
## 2 15607 16002 1 45 68 CMo 5
## 3 15608 16003 1 60 62 CMe 11
## 4 15609 16003 1 40 98 LPe 12
## 5 15618 16009 1 100 34 ANu 9
## 6 15619 16010 1 40 10 GLe 3sort(unique(records$SU_SYM90))## [1] ANu ANz ARg ARh CMd CMe CMo CMu FLd FRh GLd GLe LPd LPe LPu LVh PHh
## [18] PLe VRe
## 19 Levels: ANu ANz ARg ARh CMd CMe CMo CMu FLd FRh GLd GLe LPd LPe ... VReunique(records$MU_GLOBAL)## [1] 16002 16003 16009 16010 16013 16014 16015 16016 16017 16018 16024
## [12] 16028 16030 16031 16039 16042 16047 16049 16053 16054unique(records$SHARE)## [1] 55 45 60 40 100 20 70 30 50records[, c("ID", "MU_GLOBAL", "SHARE", "SU_SYM90", "T_SILT",
"T_GRAVEL")]## ID MU_GLOBAL SHARE SU_SYM90 T_SILT T_GRAVEL
## 1 15606 16002 55 CMd 37 0
## 2 15607 16002 45 CMo 23 3
## 3 15608 16003 60 CMe 13 1
## 4 15609 16003 40 LPe 8 19
## 5 15618 16009 100 ANu 34 0
## 6 15619 16010 40 GLe 33 0
## 7 15620 16010 40 CMe 55 0
## 8 15621 16010 20 FLd 41 3
## 9 15624 16013 60 CMo 24 0
## 10 15625 16013 40 CMu 22 0
## 11 15626 16014 70 ANu 46 0
## 12 15627 16014 30 CMd 28 0
## 13 15628 16015 70 CMu 39 1
## 14 15629 16015 30 LPd 29 9
## 15 15630 16016 70 ANu 34 0
## 16 15631 16016 30 ANz 7 0
## 17 15632 16017 60 LPu 29 1
## 18 15633 16017 40 ANu 26 0
## 19 15634 16018 60 CMu 29 0
## 20 15635 16018 40 LPd 24 22
## 21 15649 16024 70 ARg 4 0
## 22 15650 16024 30 GLd 23 0
## 23 15658 16028 60 CMd 9 0
## 24 15659 16028 40 FRh 14 4
## 25 15662 16030 50 LVh 12 0
## 26 15663 16030 50 ARh 8 0
## 27 15664 16031 55 CMe 21 4
## 28 15665 16031 45 PHh 26 1
## 29 15680 16039 40 ANu 62 0
## 30 15681 16039 40 PLe 43 0
## 31 15682 16039 20 GLd 35 0
## 32 15686 16042 55 CMe 21 4
## 33 15687 16042 45 GLd 34 0
## 34 15696 16047 60 LPd 24 22
## 35 15697 16047 40 CMu 29 0
## 36 15700 16049 100 LPu 29 1
## 37 15707 16053 50 PHh 42 0
## 38 15708 16053 30 VRe 29 0
## 39 15709 16053 20 GLd 22 0
## 40 15710 16054 70 CMe 29 0
## 41 15711 16054 30 CMd 33 0Observe que, en la región, las unidades de suelos tienen múltiples componentes.
Muchos de los atributos son datos de tipo “factor” en R, aunque en la base de datos relacional ellos eran de tipo “entero” o “caracter”. Es necesario informarle a R de eso.
Tarea 16: Convertir los atributos a datos tipo factor en R.
for (i in names(records)[c(2:5, 8:15, 17:19, 28, 45)]) {
eval(parse(text = paste("records$", i, " <- as.factor(records$",
i, ")", sep = "")))
}También se podrían asignar los nombres de los niveles de los factores usando las tablas de referencia de los metadatos (esta tarea no está implementada aún).
Tarea 17: Remover los atributos de la tabla de la región que no tengan datos.
Algunos atributos no están definidos dentro de la región. Por ejemplo, el atributo MU_SOURCE2 no se usa. Esto se puede verificar usando la función all aplicada a un vector lógico creado usando la función is.na y el operador lógico ! (“not”):
ix <- which(names(records) == "MU_SOURCE2")
all(is.na(records[, ix]))## [1] TRUELuego de encontrar esos atributos, ellos se pueden remover del dataframe, lo cual ayuda a simplificar la tabla:
df <- records
for (i in 1:length(names(records))) if (all(is.na(records[, i]))) df <- df[-i]
dim(records)## [1] 41 59dim(df)## [1] 41 49records <- df
rm(df, ix, i)Ahora tenemos una tabla únicamente con las unidades de nuestra región y únicamente con los atributos que están definidos. Esta tabla es un archivo plano que puede ser exportada para usar como hoja de cálculo o para ser importada en un programa que maneje bases de datos.
Tarea 18: Exportar la tabla de unidades de suelos como un archivo de valores separados por comas (CSV).
Para ello, podemos usar la función write.csv:
write.csv(records, file = "./HWSD_Colombia.csv")También podemos escribir directamente los archivos Excel usando la función write.xls provista por el paquete dataframes2xls. Esto tiene la ventaja que escribe correctamente los factores de R como variables de tipo caracter, no como variables de tipo entero.
require(dataframes2xls)## Loading required package: dataframes2xlswrite.xls(records, file = "./HWSD_Colombia.xls")Podemos ver los nombres de las unidades de suelo haciendo otra unión entre tablas. Para ello, podemos repetir la consulta anterior y guardar los resultados en una tabla nueva, de nombre TIEMPO. LUego, podemos usar esta tabla para realizar una unión nueva y así obtener los códigos de las unidades de suelo, los símbolos y los nombres:
sqlQuery(tcon, "SELECT * INTO TIEMPO FROM HWSD_DATA AS T
INNER JOIN NVENTANA AS U ON T.MU_GLOBAL = U.SMU_ID
ORDER BY SU_SYM90")## character(0)sqlColumns(tcon, "TIEMPO")$COLUMN_NAME## [1] "ID" "MU_GLOBAL" "MU_SOURCE1"
## [4] "MU_SOURCE2" "ISSOIL" "SHARE"
## [7] "SEQ" "SU_SYM74" "SU_CODE74"
## [10] "SU_SYM85" "SU_CODE85" "SU_SYM90"
## [13] "SU_CODE90" "T_TEXTURE" "DRAINAGE"
## [16] "REF_DEPTH" "AWC_CLASS" "PHASE1"
## [19] "PHASE2" "ROOTS" "IL"
## [22] "SWR" "ADD_PROP" "T_GRAVEL"
## [25] "T_SAND" "T_SILT" "T_CLAY"
## [28] "T_USDA_TEX_CLASS" "T_REF_BULK_DENSITY" "T_BULK_DENSITY"
## [31] "T_OC" "T_PH_H2O" "T_CEC_CLAY"
## [34] "T_CEC_SOIL" "T_BS" "T_TEB"
## [37] "T_CACO3" "T_CASO4" "T_ESP"
## [40] "T_ECE" "S_GRAVEL" "S_SAND"
## [43] "S_SILT" "S_CLAY" "S_USDA_TEX_CLASS"
## [46] "S_REF_BULK_DENSITY" "S_BULK_DENSITY" "S_OC"
## [49] "S_PH_H2O" "S_CEC_CLAY" "S_CEC_SOIL"
## [52] "S_BS" "S_TEB" "S_CACO3"
## [55] "S_CASO4" "S_ESP" "S_ECE"
## [58] "rownames" "SMU_ID"sqlColumns(tcon, "D_SYMBOL90")## TABLE_CAT
## 1 C:\\\\Users\\\\Iván Lizarazo\\Documents\\\\UN\\\\datos\\\\hwsd\\HWSD
## 2 C:\\\\Users\\\\Iván Lizarazo\\Documents\\\\UN\\\\datos\\\\hwsd\\HWSD
## 3 C:\\\\Users\\\\Iván Lizarazo\\Documents\\\\UN\\\\datos\\\\hwsd\\HWSD
## TABLE_SCHEM TABLE_NAME COLUMN_NAME DATA_TYPE TYPE_NAME COLUMN_SIZE
## 1 <NA> D_SYMBOL90 CODE 5 SMALLINT 5
## 2 <NA> D_SYMBOL90 VALUE 12 VARCHAR 50
## 3 <NA> D_SYMBOL90 SYMBOL 12 VARCHAR 4
## BUFFER_LENGTH DECIMAL_DIGITS NUM_PREC_RADIX NULLABLE REMARKS COLUMN_DEF
## 1 2 0 10 1 <NA> <NA>
## 2 100 NA NA 1 <NA> <NA>
## 3 8 NA NA 1 <NA> <NA>
## SQL_DATA_TYPE SQL_DATETIME_SUB CHAR_OCTET_LENGTH ORDINAL_POSITION
## 1 5 NA NA 1
## 2 12 NA 100 2
## 3 12 NA 8 3
## IS_NULLABLE ORDINAL
## 1 YES 1
## 2 YES 2
## 3 YES 3window.fao90 <- sqlQuery(tcon, "SELECT CODE, VALUE, SYMBOL FROM D_SYMBOL90 AS U
INNER JOIN TIEMPO AS T
ON T.SU_CODE90 = U.CODE")
head(window.fao90)## CODE VALUE SYMBOL
## 1 34 Umbric Andosols ANu
## 2 34 Umbric Andosols ANu
## 3 34 Umbric Andosols ANu
## 4 34 Umbric Andosols ANu
## 5 34 Umbric Andosols ANu
## 6 35 Vitric Andosols ANzsqlDrop(tcon, "TIEMPO")4. Mapas de atributos raster
El paquete raster no es apropiado para trabajar con bases de datos de atributos enlazadas con mapas. Para ello, es preferible utilizar el paquete sp.
Tarea 19: Convertir la ventana HWSD en un SpatialGridDataFrame, y agregar los atributos de la base de datos.
La función match encuentra la posición de un valor determinado en la tabla de referencia. Aquí se va a usar esa función para encontrar el valor SMU_ID del raster convertido que coincide con el registro en el dataframe de atributos. Luego, usamos ese índice para extraer el registro apropiado para cada celda, y para agregarlo al dataframe. Enseguida, desplegamos dos mapas de atributos: uno categórico y otro continuo.
hwsd.crop.sp <- as(hwsd.crop, "SpatialGridDataFrame")
str(hwsd.crop.sp@data)## 'data.frame': 129600 obs. of 1 variable:
## $ hwsd: int 16013 16013 16013 16013 16013 16013 16013 16013 16013 16013 ...m <- match(hwsd.crop.sp@data$hwsd, records$MU_GLOBAL)
str(m)## int [1:129600] 9 9 9 9 9 9 9 9 9 9 ...str(hwsd.crop.sp@data)## 'data.frame': 129600 obs. of 1 variable:
## $ hwsd: int 16013 16013 16013 16013 16013 16013 16013 16013 16013 16013 ...hwsd.crop.sp@data <- records[m, ]
str(hwsd.crop.sp@data)## 'data.frame': 129600 obs. of 49 variables:
## $ ID : int 15624 15624 15624 15624 15624 15624 15624 15624 15624 15624 ...
## $ MU_GLOBAL : Factor w/ 20 levels "16002","16003",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ MU_SOURCE1 : Factor w/ 20 levels "CO12","CO13",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ ISSOIL : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## $ SHARE : num 60 60 60 60 60 60 60 60 60 60 ...
## $ SEQ : int 1 1 1 1 1 1 1 1 1 1 ...
## $ SU_SYM74 : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
## $ SU_SYM85 : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
## $ SU_SYM90 : Factor w/ 19 levels "ANu","ANz","ARg",..: 7 7 7 7 7 7 7 7 7 7 ...
## $ T_TEXTURE : Factor w/ 3 levels "1","2","3": 3 3 3 3 3 3 3 3 3 3 ...
## $ REF_DEPTH : int 100 100 100 100 100 100 100 100 100 100 ...
## $ AWC_CLASS : Factor w/ 4 levels "1","2","3","5": 1 1 1 1 1 1 1 1 1 1 ...
## $ PHASE1 : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
## $ PHASE2 : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
## $ ROOTS : logi NA NA NA NA NA NA ...
## $ IL : logi NA NA NA NA NA NA ...
## $ SWR : logi NA NA NA NA NA NA ...
## $ T_GRAVEL : int 0 0 0 0 0 0 0 0 0 0 ...
## $ T_SILT : int 24 24 24 24 24 24 24 24 24 24 ...
## $ T_USDA_TEX_CLASS : Factor w/ 9 levels "1","3","4","5",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ T_BULK_DENSITY : num 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 ...
## $ T_PH_H2O : num 5 5 5 5 5 5 5 5 5 5 ...
## $ T_CEC_CLAY : num 29 29 29 29 29 29 29 29 29 29 ...
## $ T_CEC_SOIL : num 23 23 23 23 23 23 23 23 23 23 ...
## $ T_BS : num 45 45 45 45 45 45 45 45 45 45 ...
## $ T_TEB : num 10.4 10.4 10.4 10.4 10.4 10.4 10.4 10.4 10.4 10.4 ...
## $ T_CACO3 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ T_CASO4 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ T_ESP : num 2 2 2 2 2 2 2 2 2 2 ...
## $ T_ECE : num 0 0 0 0 0 0 0 0 0 0 ...
## $ S_GRAVEL : int 0 0 0 0 0 0 0 0 0 0 ...
## $ S_SAND : int 20 20 20 20 20 20 20 20 20 20 ...
## $ S_SILT : int 26 26 26 26 26 26 26 26 26 26 ...
## $ S_CLAY : int 54 54 54 54 54 54 54 54 54 54 ...
## $ S_USDA_TEX_CLASS : Factor w/ 11 levels "1","2","3","4",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ S_REF_BULK_DENSITY: num 1.22 1.22 1.22 1.22 1.22 1.22 1.22 1.22 1.22 1.22 ...
## $ S_BULK_DENSITY : num 1.19 1.19 1.19 1.19 1.19 1.19 1.19 1.19 1.19 1.19 ...
## $ S_OC : num 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 ...
## $ S_PH_H2O : num 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 ...
## $ S_CEC_CLAY : num 23 23 23 23 23 23 23 23 23 23 ...
## $ S_CEC_SOIL : num 15 15 15 15 15 15 15 15 15 15 ...
## $ S_BS : num 50 50 50 50 50 50 50 50 50 50 ...
## $ S_TEB : num 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 ...
## $ S_CACO3 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ S_CASO4 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ S_ESP : num 2 2 2 2 2 2 2 2 2 2 ...
## $ S_ECE : num 0 0 0 0 0 0 0 0 0 0 ...
## $ rownames : int 5 5 5 5 5 5 5 5 5 5 ...
## $ SMU_ID : int 16013 16013 16013 16013 16013 16013 16013 16013 16013 16013 ...rm(m)Tarea 20: Desplegar un mapa de tipos de suelo FAO 1990 y un mapa de las clases de textura en la capa superior del suelo.
Para esto podemos utilizar el método spplot, especificando la variable que se va a desplegar usando el argumento zcol:
spplot(hwsd.crop.sp, zcol = "SU_SYM90", col.regions = topo.colors(length(levels(hwsd.crop.sp$SU_SYM90))),
main = "Códigos de tipo de suelo FAO 1990", scales = list(draw = TRUE))
spplot(hwsd.crop.sp, zcol = "T_USDA_TEX_CLASS", col.regions = bpy.colors(64),
main = "Clases de textura USDA en la capa
superior del suelo",
scales = list(draw = TRUE))
5. Datos en formato vectorial
Aunque los datos HWSD están en formato raster, ellos fueron creados a partir de datos vectoriales (polígonos). Los datos vectoriales usan menos espacio de almacenamiento y usualmente son más atractivos visualmente. Los modeladores prefieren, en general, daots raster. Otros usuarios, por el contrario, prefieren datos vectoriales.
5.1 Conversión raster a vector
Tarea 21: Convertir los datos raster a una estructura vectorial, cada polígono debe tener como etiqueta el código de las celdas contiguas que lo conforman.
El paquete raster tiene la función rasterToPolygons que permite realizar esta conversión. Para usar esa función se requiere el paquete rgeos que permite disolver los límtes entre polígonos que tengan el mismo código.
require(rgeos)## Loading required package: rgeos## rgeos version: 0.3-23, (SVN revision 546)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-4
## Polygon checking: TRUEsystem.time(hwsd.crop.poly <- rasterToPolygons(hwsd.crop, n = 4,
na.rm = TRUE, dissolve = TRUE))## user system elapsed
## 60.60 1.26 61.97class(hwsd.crop.poly)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"str(hwsd.crop.poly@data)## 'data.frame': 20 obs. of 1 variable:
## $ hwsd: num 16002 16003 16009 16010 16013 ...spplot(hwsd.crop.poly, col.regions = terrain.colors(64), main = "Código de unidades de suelos",
scales = list(draw = TRUE))
Existen únicamente 20 unidades de suelos (conjuntos de polígonos con el mismo código), a diferencia de 129600 celdas raster. Esto significa unos ahorros de memoria y de tiempo de procesamiento muy grandes.
Los datos vectoriales usados por el paquete sp no están proyectados de la misma manera que los datos raster. No se requiere un remuestreo de todos los datos, solamente una re-proyección de los límites. Para ello se puede utilizar la función spTransform del paquete rgdal. Esto requiere un sistema de referencia de coordenadas (CRS) objetivo, el cual se almacena en proj4string, la cadena de especificación CRS en formato PROJ.4, en todos los objetos sp. Dicha información CRS fue almacenada en la variable proj4string.3116 al realizar la Tarea No. 6.
hwsd.crop.poly.tm <- spTransform(hwsd.crop.poly, CRS(proj4string.3116))5.2 Mapas con atributos de polígonos
Hasta ahora los datos vectoriales tienen únicamente el código de las unidades de suelo. Ahora necesitamos agregar los atributos al dataframe, y desplegar el tipo de suelo y la clase de textura en la capa superior del suelo.
El ajuste de atributos y códigos es similar al realizado en una sección previa. El nombre del atributo en los datos vectoriales es hwsd, y luego se debe comparar con el atributo MU_GLOBAL en la tabla de datos:
m <- match(hwsd.crop.poly.tm$hwsd, records$MU_GLOBAL)
str(m)## int [1:20] 1 3 5 6 9 11 13 15 17 19 ...hwsd.crop.poly.tm@data <- records[m, ]
spplot(hwsd.crop.poly.tm, zcol = "SU_SYM90", col.regions = topo.colors(length(levels(hwsd.crop.sp$SU_SYM90))),
main = "HWSD Código de unidades de suelos FAO 1990", scales = list(draw = TRUE))
spplot(hwsd.crop.poly.tm, zcol = "T_USDA_TEX_CLASS", col.regions = bpy.colors(64),
main = "Clase de textura USDA en la capa superior", scales = list(draw = TRUE))
6. Unidades de suelos con múltiples componentes
Como se observó en las secciones previas, las unidades de suelos de la zona de ejemplo tienen múltiples componentes. Revisemos nuevamente los registros obtenidos
Tarea 23: Desplegar nuevamente los componentes de las unidades de suelos.
Revisemos los múltiples componentes de las unidades de suelos en la región:
records[, c("ID", "MU_GLOBAL", "SHARE", "SU_SYM90", "T_SILT",
"T_GRAVEL")]## ID MU_GLOBAL SHARE SU_SYM90 T_SILT T_GRAVEL
## 1 15606 16002 55 CMd 37 0
## 2 15607 16002 45 CMo 23 3
## 3 15608 16003 60 CMe 13 1
## 4 15609 16003 40 LPe 8 19
## 5 15618 16009 100 ANu 34 0
## 6 15619 16010 40 GLe 33 0
## 7 15620 16010 40 CMe 55 0
## 8 15621 16010 20 FLd 41 3
## 9 15624 16013 60 CMo 24 0
## 10 15625 16013 40 CMu 22 0
## 11 15626 16014 70 ANu 46 0
## 12 15627 16014 30 CMd 28 0
## 13 15628 16015 70 CMu 39 1
## 14 15629 16015 30 LPd 29 9
## 15 15630 16016 70 ANu 34 0
## 16 15631 16016 30 ANz 7 0
## 17 15632 16017 60 LPu 29 1
## 18 15633 16017 40 ANu 26 0
## 19 15634 16018 60 CMu 29 0
## 20 15635 16018 40 LPd 24 22
## 21 15649 16024 70 ARg 4 0
## 22 15650 16024 30 GLd 23 0
## 23 15658 16028 60 CMd 9 0
## 24 15659 16028 40 FRh 14 4
## 25 15662 16030 50 LVh 12 0
## 26 15663 16030 50 ARh 8 0
## 27 15664 16031 55 CMe 21 4
## 28 15665 16031 45 PHh 26 1
## 29 15680 16039 40 ANu 62 0
## 30 15681 16039 40 PLe 43 0
## 31 15682 16039 20 GLd 35 0
## 32 15686 16042 55 CMe 21 4
## 33 15687 16042 45 GLd 34 0
## 34 15696 16047 60 LPd 24 22
## 35 15697 16047 40 CMu 29 0
## 36 15700 16049 100 LPu 29 1
## 37 15707 16053 50 PHh 42 0
## 38 15708 16053 30 VRe 29 0
## 39 15709 16053 20 GLd 22 0
## 40 15710 16054 70 CMe 29 0
## 41 15711 16054 30 CMd 33 0Por ejemplo, la unidad 16010 (cuyos registros son 15619, 15620 y 15621) tiene tres componentes, con porcentajes de 40%, 40% y 20%. Los tres componentes tienen la siguiente concentración de limos: 33%, 55% y 41%. Por otra parte, la unidad 16009 solamente tiene un componente y su concentración de limo es 34%.
Las unidades cartográficos con más de un componente crean problemas para elaborar mapas de atributos en formato raster. Esto significa que la técnica explicada en la sección 4 debe ser modificada pues más de un registro puede cumplir la unión de tablas. Existen varias soluciones a este problema. Al usar la función match se encuentra el primer componente de la lista, de manera tal que todos los atributos de una unión simple corresponden únicamente a ese primer componente. Por ejemplo, el contenido de limos en la parte superior del suelo:
hwsd.crop.sp <- as(hwsd.crop, "SpatialGridDataFrame")
m <- match(hwsd.crop.sp@data$hwsd, records$MU_GLOBAL)
hwsd.crop.sp@data <- records[m, ]
summary(hwsd.crop.sp@data$T_SILT)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 24.00 29.00 31.49 39.00 62.00Al comparar este resultado con la lista completa de componentes de las unidades de suelo, se puede observar que solamente se lista el primer componente, por ejemplo, el contenido de limo es solamente para ese componente. En consecuencia, un mapa muestra únicamente el atributo del primer componente de la lista: A map shows the attribute of the first-listed component:
spplot(hwsd.crop.sp, zcol = "T_SILT", col.regions = bpy.colors(64),
main = "Proporcion de limo en la capa superior,
%, suelo dominante",
scales = list(draw = TRUE))
Tarea 24: Promediar el contenido de limo
Existen varias opciones para calcular un valor para cada celda:
Aceptar el valor obtenido mediante la función match, es decir, el valor correspondiente al primer componente de la lista. Si los valores están listados en orden descendente, este valor debería corresponder al valor dominante. Sin embargo, este ordenamiento no se puede garantizar completamente.
Seleccionar un valor único que puede ser, el valor más alto, el más bajo, la mediana, de acuerdo con la aplicación. Por ejemplo, en un mapa de riesgo de contaminación de agua subterránea se podría seleccionar el valor más alto del contenido de arena.
Calcular un valor promedio asignados pesos en función de la proporción. Para ello, se puede usar un operador SQL “aggregrate”: MAX, MIN, AVG. Para ello se requiere una cláusula adicional en la declaración SQL, es decir SQL GROUP BY, con el objeto de agregar primero los registros que están relacionados y, luego, aplicar la función. Estas funciones pueden también usar el modificador AS para renombrar el atributo resultante.
Por ejemplo, para promediar el contenido de limo:
TMP1 <- sqlQuery(tcon, "select T.ID, T.MU_GLOBAL, U.SMU_ID, T.SU_SYM90,
T.T_SILT from HWSD_DATA as T inner
join NVENTANA as U on T.MU_GLOBAL=U.SMU_ID
order by T.MU_GLOBAL")
print(TMP1)## ID MU_GLOBAL SMU_ID SU_SYM90 T_SILT
## 1 15607 16002 16002 CMo 23
## 2 15606 16002 16002 CMd 37
## 3 15608 16003 16003 CMe 13
## 4 15609 16003 16003 LPe 8
## 5 15618 16009 16009 ANu 34
## 6 15619 16010 16010 GLe 33
## 7 15620 16010 16010 CMe 55
## 8 15621 16010 16010 FLd 41
## 9 15625 16013 16013 CMu 22
## 10 15624 16013 16013 CMo 24
## 11 15626 16014 16014 ANu 46
## 12 15627 16014 16014 CMd 28
## 13 15628 16015 16015 CMu 39
## 14 15629 16015 16015 LPd 29
## 15 15630 16016 16016 ANu 34
## 16 15631 16016 16016 ANz 7
## 17 15632 16017 16017 LPu 29
## 18 15633 16017 16017 ANu 26
## 19 15635 16018 16018 LPd 24
## 20 15634 16018 16018 CMu 29
## 21 15649 16024 16024 ARg 4
## 22 15650 16024 16024 GLd 23
## 23 15658 16028 16028 CMd 9
## 24 15659 16028 16028 FRh 14
## 25 15662 16030 16030 LVh 12
## 26 15663 16030 16030 ARh 8
## 27 15664 16031 16031 CMe 21
## 28 15665 16031 16031 PHh 26
## 29 15681 16039 16039 PLe 43
## 30 15682 16039 16039 GLd 35
## 31 15680 16039 16039 ANu 62
## 32 15687 16042 16042 GLd 34
## 33 15686 16042 16042 CMe 21
## 34 15696 16047 16047 LPd 24
## 35 15697 16047 16047 CMu 29
## 36 15700 16049 16049 LPu 29
## 37 15707 16053 16053 PHh 42
## 38 15708 16053 16053 VRe 29
## 39 15709 16053 16053 GLd 22
## 40 15711 16054 16054 CMd 33
## 41 15710 16054 16054 CMe 29s_average <- aggregate(TMP1[, 5], list(TMP1$MU_GLOBAL), mean)
print(s_average)## Group.1 x
## 1 16002 30.00000
## 2 16003 10.50000
## 3 16009 34.00000
## 4 16010 43.00000
## 5 16013 23.00000
## 6 16014 37.00000
## 7 16015 34.00000
## 8 16016 20.50000
## 9 16017 27.50000
## 10 16018 26.50000
## 11 16024 13.50000
## 12 16028 11.50000
## 13 16030 10.00000
## 14 16031 23.50000
## 15 16039 46.66667
## 16 16042 27.50000
## 17 16047 26.50000
## 18 16049 29.00000
## 19 16053 31.00000
## 20 16054 31.00000Compare el resultado anterior con la tabla de unidades de suelos. Observe que las unidades de suelo que tienen múltiples componentes se han agregado. Por ejemplo, la unidad 16010, que tenía tres registros, muestra ahora el promedio del contenido de limo. Este resultado no es exactamente lo que buscamos debido a que:
no siempre las proporciones de los componentes son iguales;
si uno de los componentes no tiene limo, el promedio debería calcularse usando este cero implícito.
Para obtener un promedio ponderado correcto, se requiere extraer las proporciones y el peso del contenido de limo. Por ahora, esta actividad se deja como ejercicio para el lector.
7. Limpieza
Tarea 25 : Remover las tablas temporales y desconectar la base de datos
# sqlDrop(tcon, 'WINDOW_CROP')
sqlDrop(tcon, "NVENTANA")
odbcClose(tcon)Referencias
[1] Roger S. Bivand, Edzer J. Pebesma, and V. Gómez-Rubio. Applied Spatial Data Analysis with R. Springer, 2008. ISBN 978-0-387-78170- 9. URL http://www.asdar-book.org/.
[2] Robert J. Hijmans and Jacob van Etten. raster: Geographic analysis and modeling with raster data, 2012. URL http://CRAN.R-project. org/package=raster. R package version 2.0-12.
[3] IIASA; FAO; ISRIC; ISS-CAS; JRC. Harmonized World Soil Database (version 1.2). FAO and IIASA, Feb 2012. URL http://webarchive. iiasa.ac.at/Research/LUC/External-World-soil-database/ HWSD_Documentation.pdf.
[4] Timothy H. Keitt, Roger Bivand, Edzer Pebesma, and Barry Rowlingson. rgdal: Bindings for the Geospatial Data Abstraction Library, 2012. URL http://CRAN.R-project.org/package=rgdal. R package version 0.7-20.
[5] osgeo.org. PROJ.4. URL https://trac.osgeo.org/proj/.
[6] Edzer J. Pebesma and Roger S. Bivand. Classes and methods for spatial data in R. R News, 5(2):9-13, 2005.
[7] R Development Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2012. URL http://www.R-project.org/. ISBN 3-900051- 07-0.
[8] R Development Core Team. R Data Import/Export. The R Foundation for Statistical Computing, version 2.15.2 (2012-10-26) edition, 2012. ISBN ISBN 3-900051-10-0. URL http://cran.r-project.org/doc/ manuals/R-data.pdf.
[9] Yihui Xie. knitr: Elegant, flexible and fast dynamic report generation with R, 2011. URL http://yihui.name/knitr/.
[10] Brian Ripley and Michael Lapsley. Package ‘RODBC’, 2017. URL https://cran.r-project.org/web/packages/RODBC/RODBC.pdf