Conversión de Archivos KML a Shapefile mediante procesamiento en Python y R

En este tutorial explicare el proceso de importar archivos KML a shapefile mediante dos códigos, el primero escrito en el lenguaje de programación Python, y el segundo en el lenguaje R; en el primer proceso se exportan cada archivo Kml a una geodatabase, y luego en R importamos los archivos de la geodatadase y los exportamos como un archivo shapefile. Para llevar a cabo el presente tutorial es necesario tener instalado cualquier versión de ArcMAP superior a 10, y RStudio versión 3 o superior.

La ventaja que nos brinda hacer el proceso mediante línea de código es que si tenemos, por ejemplo, 100 archivos KML no sería necesario hacer el proceso 100 veces con la función ‘Import KML to Layer’ de ArcGIS, sino de una manera mucho más eficiente con un ciclo en ambos lenguajes de programación.

Los archivos para realizar el presente tutorial están disponibles en este enlace web. Estos archivos corresponden a polígonos dibujados aleatoriamente en Google Earth, no corresponden a ningún tipo de cobertura, uso de la tierra, o elemento geográfico; fueron realizados en Google Earth solamente para fines de hacer este tutorial.

Aquí los pasos a seguir:

Sección Python

1. En un compilador de texto, escribimos las siguientes líneas de código, el compilador que yo utilizare será notepad (clic aquí para descargar). Las primeras líneas son las siguientes:

import arcpy, os
arcpy.env.workspace = r'path\_kmz'
outLocation = r'path\_gdb'
MasterGDB = 'AllKMLLayers.gdb'
MasterGDBLocation = os.path.join(outLocation, MasterGDB)

En la segunda línea indicamos cuál será nuestro espacio de trabajo, en este espacio de trabajo es donde se encuentra el listado de archivos KMZ. En la tercer linea va la ruta donde quedaran guardadas las geodatabases. Y por último, en la cuarta línea escribimos el nombre de cómo queremos que llame la geodatabase donde guardaremos los layers.
2.      Ahora escribiremos estas líneas, aquí lo que haremos es hacer un for para leer los archivos .KM que contiene la carpeta de nuestro espacio de trabajo, y establecemos el ambiente de trabajo donde se encontraran nuestros archivos de salida.
print 'Master'

for kmz in arcpy.ListFiles('*.KM*'):
print ("CONVERTING: {0}".format(os.path.join(arcpy.env.workspace, kmz)))
arcpy.KMLToLayer_conversion(kmz, outLocation)

arcpy.env.workspace = outLocation
wks = arcpy.ListWorkspaces('*', 'FileGDB')
wks.remove(MasterGDBLocation)
1.   3.   El siguiente paso será escribir el último grupo de líneas, en donde se ejecuta la conversión de los archivos kmz a geodatabase mediante un for, este permite iterar en caso dado que se tengan muchos archivos a convertir

for fdgb in wks:

arcpy.env.workspace = fdgb
featureClasses = arcpy.ListFeatureClasses('*', '', 'Placemarks')

for fc in featureClasses:
print ("COPYING: {0} FROM: {1}".format(fc, fgdb))
fcCopy = os.path.join(fgdb, 'Placemarks', fc)
arcpy.FeatureClassToFeatureClass_conversion(fcCopy, MasterGDBLocation, fgdb[fgdb.rfind(os.sep)+1:-4])


# Clean up
del kmz, wks, fc, featureClasses, fgdb

Estas líneas de código fueron escritas a partir de la ayuda que ofrece ArcMap (sí desea visitar la página web aquí el enlace).
Sección R
Ahora pasaremos a abrir RStudio y escribiremos las siguientes líneas.
1.   1. El primer pasos será escribir, de igual manera como en Python, la carpeta de trabajo en donde se encuentran las geodatabases creadas mediante el punto anterior. Antes de ello es importante cargar las librerías a utilizar, en caso dado que usted no tenga las librerías las puede instalar mediante el comando install.packages (‘nameLibrary’).
require(raster)
require(rgdal)
require(plotKML)
require(maptools)
require(XML)

files <- list.files('_data/_gdb', full.names = T, pattern = '.gdb')
list <- lapply(files, FUN = ogrListLayers)

En la línea del lapply aplicamos la función ogrListLayers a las geodatabases, esto con la finalidad de conocer cuál es el nombre de cada shapefile dentro de la geodatabase.
2. Luego creamos un objetos tipo lista que es donde guardaremos los archivos con los shapefiles
poly  <- list()
1.  
      3. Seguido de ello pasamos a escribir un ciclo for, mediante el cuál leemos los archivos shapefiles dentro de la geodatabase, y escribimos estos archivos en una carpeta distinta a la que contiene las geodatabases. Dentro del for hay un condicional, este sirve para conocer si el archivo que vamos a escribir ya existe dentro de nuestro espacio de trabajo, en caso dado que ya exista el archivo no lo escribe de nuevo. 

for(i in 1:length(files)){

poly[[i]] <- readOGR(dsn = files[[i]], layer = as.character(list[[i]][1]))

if(!file.exists(paste0('/_data/_shp/', gsub(' ', '', as.character(poly[[i]]@data[1,1])), '.shp'))){

writeOGR(obj = poly[[i]], dsn = '_data/_shp', layer = gsub(' ', '', as.character(poly[[i]]@data[1,1])), driver = 'ESRI Shapefile')
print('Written Shapefiles')

}

}

Los códigos los puedes encontrar, además, en mi cuenta de Github, disponible en este enlace.

Realizar extracción por máscara en R haciendo uso de paralelización para optimizar el tiempo en el procesamiento

En este tutorial realizaremos una extracción por máscara de 48 archivos raster, que representan la precipitación, temperatura máxima, media y mínima, para 3 distintos departamentos de Colombia; la idea final es obtener los 48 raster para cada departamento, es decir que se obtendrían al final 144 raster. Para ello haremos uso de paralelización, lo que permite ejecutar el mismo proceso para los tres distintos departamentos en núcleos distintos, es decir, a la misma vez que se cortan los datos para Antioquia se están cortando los datos para Cundinamarca y Valle del Cauca, pues se hará uso de tres núcleos (si su computador tiene menos de cuatro núcleos no podrá replicar este ejercicio tal cual).

Los tres departamentos a los cuales les haremos la extracción por máscara son Antioquia, Valle del Cauca y Cundinamarca, para ello tendremos el shapefile de división político administrativa de Colombia (Fuente: SIGOT); además de ello, los raster climáticos para toda Colombia. El ejercicio incluye la descarga de los datos para que replique el ejercicio (URL).

1. Cargamos las librerías necesarias, en caso dado que usted no tenga la librería esta se descargara automáticamente.

if(!require(raster)){install.packages('raster'); library(raster)} else {library(raster)}
if(!require(rgdal)){install.packages('rgdal'); library(rgdal)} else {library(rgdal)}
if(!require(spdplyr)){install.packages('spdplyr'); library(spdplyr)} else {library(spdplyr)}
if(!require(tidyverse)){install.packages('tidyverse'); library(tidyverse)} else {library(tidyverse)}
if(!require(foreach)){install.packages('foreach'); library(foreach)} else {library(foreach)}
if(!require(parallel)){install.packages('parallel'); library(parallel)} else {library(parallel)}
if(!require(doSNOW)){install.packages('doSNOW'); library(doSNOW)} else {library(doSNOW)}
if(!require(gridExtra)){install.packages('gridExtra'); library(gridExtra)} else {library(gridExtra)}
if(!require(rasterVis)){install.packages('rasterVis'); library(rasterVis)} else {library(rasterVis)}

2. Descargamos y declaramos los archivos a utilizar en este tutorial, tanto los archivos raster como el shapefile de departamentos de Colombia

url    <- 'https://zdl9fa.bn1303.livefilestore.com/y4mV948RDpIKmA015QmtTAON46y0O8vW0kF7B3JQ6xeP006WJDddpGmwQ1wpRS9G3f36LGwVNy-EbwxqA9HPNUX6msaosaUNe3RJf8WynIjTbMohC490Mer7l5lrpEMTdVPR1-XWmPHv9CUaPxLFnqqu1qLIo5a36EjDeR4Ry0QmCsvt4G9TXyoLVoMkWt54VOUpAH4VIsLSuS5yk44wVI8Iw/_data.zip'
mydir  <- 'D:/parallel'
if(!file.exists(mydir)){dir.create(mydir)}
temp <- tempfile(tmpdir = mydir, fileext = '.zip')
download.file(url, temp)
unzip(temp, exdir = mydir)
unlink(temp) # Eliminar el archivo .zip

adm <- shapefile(paste0(mydir, '/_data/_shp/LimiteDptal_geo_ok.shp'))
lyrs <- list.files(paste0(mydir, '/_data/_raster/_co'), full.names = T, pattern = '.asc$') %>%
lapply(FUN = raster) %>%
stack()

3.  Del shape de Colombia extraemos tres objetos distintos que sean cada departamento, y luego los guardamos en un objeto tipo lista. 

ant   <- filter(adm, NOMBRE_DPT == 'ANTIOQUIA')
val <- filter(adm, NOMBRE_DPT == 'VALLE DEL CAUCA')
cun <- filter(adm, NOMBRE_DPT == 'CUNDINAMARCA')
shps <- list(ant, val, cun)

4. Creamos las carpetas donde se guardaran los archivos cortados para cada departamento, en mi caso creare ‘_ant’ para Antioquía, ‘_val’ para Valle del Cauca, ‘_cun’ para Cundinamarca.

dirsOut <- paste0(path, '/', c('_ant', '_val', '_cun'))
if(!file.exists(dirsOut)){lapply(dirsOut, FUN = dir.create)}

5.  Indicamos la cantidad de núcleos a utilizar. Se recomienda utilizar un núcleo menos que el total de núcleos que tenga su computador.

nCores   <- detectCores(all.tests = FALSE, logical = TRUE)
cl <- makeCluster(nCores - 1) #Número de nucleos a utilizar
registerDoSNOW(cl)

6. Configuramos la barra de progreso para así conocer que tanto hace falta para que se termine el proceso de extracción por máscara.

length_run <- length(shps)
pb <- txtProgressBar(max = length_run, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress=progress)


7. Creación de la función de extracción por máscara. Primero usamos la función, crop y mask para cortar por departamento y luego realizamos un unstack para tener los archivos raster como lista, y poderlos así escribir mediante un for sencillo.

extractMask <- function(lyrs, shps, dirsOut){

lyrs_cut <- crop(lyrs, shps) %>%
mask(shps) %>%
unstack()

for(j in 1:length(lyrs_cut)){

writeRaster(lyrs_cut[[j]], paste0(dirsOut, '/', names(lyrs_cut[[j]]), '.asc'))

}

return(lyrs_cut)

}

8. Ejecutamos el foreach para correr el proceso paralelamente en los tres núcleos destinados para ello, se recomienda no tener abierto muchos programas en el computador para que así no se ponga lento el PC.

lyrs_cut_ok <- foreach(i = 1:length(shps), .packages = c('raster', 'rgdal', 'dplyr'), .options.snow=opts) %dopar% {

extractMask(lyrs, shps[[i]], dirsOut[i])

}

9. Creamos una función sencilla para plotear los archivos resultantes

plotGrid <- function(lyr){

gg <- gplot(lyr) +
geom_tile(aes(fill = value)) +
coord_equal() +
labs(fill = "Precipitación \n Enero (mm)") +
xlab('Lon') +
ylab('Lat')

}

plotAnt <- plotGrid(lyr = lyrs_cut_ok[[1]][[1]])
plotCun <- plotGrid(lyr = lyrs_cut_ok[[2]][[1]])
plotVal <- plotGrid(lyr = lyrs_cut_ok[[3]][[1]])

plotAll <- grid.arrange(plotAnt, plotCun, plotVal, ncol=3)

Aquí podemos observar la precipitación para el mes de enero en cada departamento:

Muchas gracias a Jeison Mesa por la ayuda brindada en la construcción de este tutorial.

Función Zonal Histogram de ArcGIS en R

En este tutorial explicare una aplicación de la función “zonal histogram” de ArcGIS en R, el objetivo, será calcular el área de las coberturas para los municipios del Valle del Cauca a partir de la capa de cobertura de la NASA y la división político administrativa del SIGOT.
Los datos a utilizar, en coordenadas geográficas WGS 1984, son los siguientes (disponibles aquí):
1.      Cobertura de la tierra reclasificada de la NASA.
2.      Municipios del Valle del Cauca
A continuación los pasos:
1.      Cargamos las librerías a utilizar, en caso que usted no las tenga instalada podría utilizar install.packages(“nameLibrary”)

# Cargamos librerias
require(raster)
require(rgdal)
require(foreign)
require(ggplot2)
require(dplyr)
2.      Cargamos los archivos, que son el raster, la tabla de la leyenda del raster, y el shapefile de municipios.

path <- "E:/R/_tutoriales/_tabulateArea/_blog/_data"
cob <- raster(paste0(path, "/_raster/cov_eluglcde.tif"))
adm1 <- shapefile(paste0(path, "/_shp/mpiosValle.shp"))
leg <- read.dbf(paste0(path, "/_raster/cov_eluglcde.tif.vat.dbf"))

3.      Ploteamos el shapefile para observar los límites de los municipios.

ggplot() + geom_polygon(data=adm1,  aes(x=long, y=lat, group=group), fill="cadetblue", color="grey") + coord_equal()
4.      Creamos la función de zonal histogram, pues esta función aún no existe como tal en alguna librería.

tabFunc                            <- function(indx, extracted, region, regname) {
dat <- as.data.frame(table(extracted[[indx]]))
dat$name <- region[[regname]][[indx]]
return(dat)
}
5.      Ejecutamos la función creada anteriormente

ext <- extract(cob, adm1, method='simple')#extract valores
tabs <- lapply(seq(ext), tabFunc, ext, adm1, "NOM_MUNICI")
df <- do.call(rbind, tabs)#unimos todas las tablas en una sola
6.      Calculamos el porcentaje de cobertura por municipio y las hectáreas, teniendo en cuenta que cada píxel tiene un área aproximada de 0.20833 hectáreas

df  %
group_by(name)%>%
mutate(Porcentaje = (Freq/sum(Freq))*100)%>%
ungroup()%>%
mutate(HA = Freq * 0.2083333)

7.      Editamos los encabezados y unimos con la tabla de leyenda para conocer cada valor de celda a que categoría corresponde ej. El valor 9 corresponde a áreas urbanas.

colnames(df)      <- c("Category", "Count", "Municipio", "Porcentaje", "Has")
df_all <- merge(df, leg, by.x = "Category", by.y = "Value")
df_all <- df_all[,c("Category", "Count.x", "Municipio", "Porcentaje", "Has", "ELU_GLC_De")]
df_all <- df_all[order(df_all$Municipio, decreasing = FALSE), ]
colnames(df_all) <- c("Category", "Count", "Municipio", "Porcentaje", "Has", "Count2", "Cobertura")
      8.       Creamos gráfica para tres municipios del Valle del Cauca, esta gráfica representa el porcentaje de cada cobertura en el municipio según el raster de cobertura de la NASA para el año 2015.

mpios <- filter(df_all, Municipio %in% c("YUMBO", "CALI", "PALMIRA"))
gg <- ggplot(mpios, aes(x = factor(Municipio), y = Porcentaje, group = Cobertura)) +
geom_bar(aes(fill = Cobertura), stat = "identity", position = "dodge") +
theme(legend.position = "bottom") + xlab("Municipio")
ggsave(plot = gg, paste0(path, "/cob_mpios_porc.png"), width = 9, height = 8, units = "in")
gg
dev.off()


9.      Guardamos la tabla como un archivo csv

write.csv(df_all, paste0(path, "/_table/zonalHistogram.csv"))

Gracias a Jeison Mesa por su colaboración en la realización de este tutorial!

Hacer gráfica en R de precipitación a partir de archivos raster por departamento


En este tutorial explicare brevemente como realizar un gráfico tipo ojiva, en el cual se presente en el eje x los meses del año & en el eje y la precipitación acumulada para cada mes, a tres distintos departamentos de Colombia, estos son: Antioquía, Cundinamarca y Valle del Cauca. Los datos de entrada fueron tomados de la base de datos Worldclim (Hijmans et al., 2005) y extraídos para Colombia a partir del shape de límite departamental del SIGOT.
Los datos utilizados en este tutorial están disponibles en este enlace, y se describen a continuación:
1) 12 raster de precipitación para todo Colombia a resolución de 30 arcos de segundo (aproximadamente 1 km2)
2) Shapefile de los 32 departamentos de Colombia.
A continuación se describen los pasos, asumiendo que ya se tiene instalado R y RStudio en su computador.
1. Instalamos las siguientes librerías (importante tener conexión a Internet)
install.packages("raster")
install.packages("rgdal")
install.packages("ggplot2")
install.packages("dplyr")
install.packages("reshape2")

require("raster")
require("rgdal")
require("ggplot2")
require("dplyr")
require("reshape2")
2.  Llamamos nuestros datos en RStudio

myproj <- CRS("+proj=longlat +datum=WGS84") #proyección a utilizar
path <- "E:/Blogs/_blogger/_graphClimate" #dirección donde está el directorio raíz de los datos
col_shp <- shapefile(paste0(path, "/_shp/LimiteDptal.shp")) # limite departamental
files <- list.files(paste0(path, "/_climate/_colombia/_asc"), full.names = T, pattern = ".asc$") #archivos raster
 3. Escogemos en nuevos objetos los departamentos de interés

valleCauca_shp <- col_shp[col_shp$NOMBRE_DPT %in% "VALLE DEL CAUCA",]
cundinamarca_shp <- col_shp[col_shp$NOMBRE_DPT %in% "CUNDINAMARCA",]
antioquia_shp <- col_shp[col_shp$NOMBRE_DPT %in% "ANTIOQUIA",]
4. Archivos raster

files_prec <- grep("prec", files, value = T) %>%
mixedsort() #seleccion de precipitación y ordenar datos
layers_prec <- lapply(files_prec, raster) #conversión de objetos a raster
5. Corte de archivos raster por departamento, se crean listas vacías donde se guardaran los raster para cada departamento.
prec_antioquia    <- list() #creación de lista vacía
prec_cundinamarca <- list() #creación de lista vacia
prec_valle <- list() #creación de lista vacia

for(i in 1:length(layers_prec)){

prec_antioquia[[i]] <- crop(layers_prec[[i]], antioquia_shp)
prec_antioquia[[i]] <- mask(prec_antioquia[[i]], antioquia_shp)

prec_cundinamarca[[i]] <- crop(layers_prec[[i]], cundinamarca_shp)
prec_cundinamarca[[i]] <- mask(prec_cundinamarca[[i]], cundinamarca_shp)

prec_valle[[i]] <- crop(layers_prec[[i]], valleCauca_shp)
prec_valle[[i]] <- mask(prec_valle[[i]], valleCauca_shp)

}
6. Calculo de promedios de precipitación para cada mes en nuestros departamentos objeto de estudio. Se crean listas vacías donde se guardaran los datos de promedios para cada raster por departamento.

prec_mean_antioquia <- list(); prec_mean_cund <- list(); prec_mean_valle <- list()
prec_df_antioquia <- NA; prec_df_cund <- NA; prec_df_valle <- NA

for(i in 1:length(prec_antioquia)){

prec_mean_antioquia[[i]] <- mean(prec_antioquia[[i]][], na.rm = T)
prec_df_antioquia <- t(as.data.frame((rbind(prec_mean_antioquia))))

prec_mean_cund[[i]] <- mean(prec_cundinamarca[[i]][], na.rm = T)
prec_df_cund <- t(as.data.frame((rbind(prec_mean_cund))))

prec_mean_valle[[i]]<- mean(prec_valle[[i]][], na.rm = T)
prec_df_valle <- t(as.data.frame((rbind(prec_mean_valle))))

}
7. Orden de datos para ggplot

df_prec <- data.frame(prec_df_antioquia, prec_df_cund, prec_df_valle)
df_prec %
mutate(Month_numeric = 1:nrow(df_prec)) %>%
mutate(prec_mean_antioquia = as.numeric(prec_mean_antioquia),
      prec_mean_cund = as.numeric(prec_mean_cund),
      prec_mean_valle  = as.numeric(prec_mean_valle))

prec_df_rename <- rename(df_prec, Antioquia = prec_mean_antioquia,
Cundinamarca = prec_mean_cund, Valle = prec_mean_valle)

prec_df_melt <- melt(prec_df_rename, id = 'Month_numeric')
prec_df_melt <- rename(prec_df_melt, Departamento = variable)
8. Creación de gráfico base

gg  <- ggplot(prec_df_melt, aes(Month_numeric, value,  fill = Departamento, colour = Departamento))+ geom_line() + xlab("Mes") + ylab("Precipitación (mm)") + ggtitle("Precipitación para cada mes") +
scale_x_continuous(expand = c(0,0), breaks = c(1:12), labels = c("Enero", "Febrero", "Marzo", "Abril", "Mayo", "Junio", "Julio", "Agosto", "Septiembre", "Octubre", "Noviembre", "Diciembre")) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(legend.position = "bottom")
9. Creación de gráfico base

gg <- gg + scale_colour_manual(labels = c("Antioquia", "Cundinamarca", "Valle del Cauca"), values = c("red", "green", "blue"))

10. Guardado de gráfico en buena resolución

ggsave(plot = gg, paste0(path, "/_png/precYear_dptos.png"), width = 9, height = 8, units = "in")
gg
El gráfico debe quedar así:

Importar y exportar archivos GPX en R

El objetivo del presente tutorial es aprender a importar un listado de archivos GPX (son los archivos nativos que se toman con el GPS) en R y exportarlos como archivos .csv, para aprender a exportar estos archivos csv a shapefile se puede leer este tutorialrealizado meses atrás.

Primero instalamos la librería “plotKML” y “data.table”, mediante el siguiente comando:

install.packages(“plotKML”)
install.packages(“data.table”)
require(plotKML)
require(data.table)
Escribimos la ruta donde está el listado de archivos gpx.
path <- "D:/Personal/_blog/_importGPX/_gpx"
listado <- list.files(path, full.names = T)
Ahora mediante el comando “lapply” aplicamos la función readGPX a todos los archivos gpx.
gpx_wayPoints_2 <- lapply(listado, FUN = readGPX, waypoints = TRUE, track = FALSE, routes = FALSE)
Ahora creamos un listado vacío para allí guardar las tablas de waypoints (puntos tomados con el GPS)
gpx_wayPoints3 <- list()
Para obtener solamente las tablas del archivo gpx, realizamos un ciclo mediante el cual iteramos en cada objeto de la lista obteniendo solamente la tabla
for (i in 1:length(listado)){
  gpx_wayPoints3[[i]] <- gpx_wayPoints_2[[i]]$waypoints
}
Mediante la siguiente línea de comando podemos observar las dimensiones de cada tabla, y así poder observar si tienen igual número de columnas
lapply(gpx_wayPointsList2, function(x){dim(x)})
Ahora mediante el comando “rbindlist” unimos los tres dataframe (tabla) dentro de uno solo
tabla <- rbindlist(gpx_wayPoints3)
dim(tabla) #podemos observer las dimensiones de la tabla
Por último guardamos la tabla en un archivo csv.
write.csv(tabla, “D:/Personal/_blog/_importGPX/tabla.csv”)
Los datos para realizar la práctica se encuentran disponibles en este enlace.

Descarga de automatizada de archivos raster en R

26 de Julio de 2016
Hace pocos días ISRIC World Soil Information [1] realizó la actualización de los datos espaciales (tipo raster) de suelo a resolución de 250 metros, anteriormente la resolución estaba en aproximadamente 1 km, estos datos provienen de distintas fuentes de datos y representan propiedades físicas, biológicas y químicas; por ejemplo, el potencial de hidrogeno (pH), la capacidad de intercambio catiónico, el contenido de carbono orgánico, el porcentaje de arena, suelo y arcilla, la densidad aparente, entre otras.
Los archivos se encuentran disponibles mediante un sitio web FPT (File Transfer Protocol) disponible en este enlace; estos datos es posible descargarlos mediante clic en cada archivo desde un navegador web (ej. Google Chrome, Opera, Mozilla Firefoz) o desde el lenguaje de programación R usando la librería “utils” (R Core Team, Versión 3.2.2.), la ventaja que ofrece R respecto al navegador web es que es posible optimizar la descarga haciendo uso de iteraciones, lo que evita el hacer uso de clic en cada archivo del sitio web [2].

install.packages(“utils”)
library(utils)
for (nomen in c(“CECSOL”,”CLYPPT”,”CRFVOL”,”ORCDRC”,”PHICKL”,”PHIHOX”, “SLTPPT”,”SNDPPT”,”ORCDRC”)) #se varia acorde al tipo de variable
{
  for (i in 1:7)
  {
    code=paste0(nomen,”_M_sl”,i,”_250m_ll.tif”)
    url=paste0(“ftp://ftp.soilgrids.org/data/recent/”,code)
    download.file(url, destfile = paste0(“E:/ISRIC/”,code), mode = “wb”)
  } #la URL es el enlace web desde donde se descargan los archivos raster y varia acorde al perfil de la variable, el “destfile” es el archivo donde se van a guardar los datos en el disco duro del pc y el nombre es el mismo del raster, el mode es wb para que sea reconocido como archivo espacial
}

En este sentido, el objetivo del presente tutorial es realizar la descarga de la mayoría de archivos raster [3] de suelos a resolución de 250 metros mediante un ciclo (for) en R, para ello necesitamos la librería “utils” y una buena conexión a internet. 
Apenas se ejecuta el código se comienza la descarga desplegándose lo presentado en la siguiente figura en caso que se esté utilizando la interfaz de RStudio.
Para tener en cuenta:
[1] ISRIC es un sistema para el mapeo automatizado del suelo basado en perfiles de suelos globales, la base de datos comprende archivos raster del mundo a 1 km y 250 m de resolución espacial producida mediante la asignación automática del suelo basada en algoritmos de aprendizaje automático. Para leer más: Hengl, T., Mendes de Jesus, J., Heuvelink, G. B.M., Ruiperez Gonzalez, M., Kilibarda, M. et al. (2016) SoilGrids250m: global gridded soil information based on Machine Learning. Earth System Science Data (ESSD) in review.
[2] Desde cualquier navegador web solamente es posible realizar dos descargas a la vez.
[3] Debido a los nombres de los archivos raster de tipo de suelo no es posible realizar la descarga optimizada mediante el ciclo, pues el nombre del archivo varía considerablemente.
Gracias @AndresAgui90 por su grata ayuda en la elaboración de este tutorial. 

Uso de la herramienta “Clip” en ArcGIS

17 de enero de 2016

Uso de la herramienta “Clip” en ArcGIS


El objetivo de este ejercicio es obtener los municipios solamente del Valle del Cauca a partir del shape de municipios para Colombia y el límite de este departamento, estos datos son descargados de la página oficial del SIGOT.
Nota: el tutorial presente fue realizado con el software ArcGIS 10.3, de igual manera la herramienta se encuentra en versiones anteriores.

Primero abrimos el shapefile de municipios y límite del Valle del Cauca en ArcGIS usando el icono resaltado en amarillo: 
Deberemos tener en el espacio de trabajo de ArcGIS algo como esto:
La herramienta de la cual se hará uso es “Clip”, esta se encuentra dentro del “ArcToolbox à Analysis Tools à Extraction à Clip”.
Se desplegara una ventana, en el archivo de entrada va a ir los municipios, en el campo “Clip Features” va el shape del departamento del Valle del Cauca, en el campo “output” va el archivo de salida (es decir, donde queremos el nuevo shape que se creara), algo importante es que ambas capas deben estar bajo el mismo sistema de coordenadas.

El archivo generado debe ser algo similar a lo siguiente:
Otra manera de obtener solamente los municipios para el Valle del Cauca es hacer una selección por atributos al shape de Municipios desde la tabla de atributos (Select by Attributes), como se muestra a continuación:
Luego damos clic derecho sobre la capa en la tabla de contenidos, luego en “Data” y “Export Data”, y guardamos el shapefile en un sitio de interés de nuestro disco duro.
Los archivos usados en el presente tutorial se encuentra disponible en este link.

Uso de la herramienta “Dissolve” en ArcGIS

10 de enero de 2016

Uso de la herramienta “Dissolve” en ArcGIS

El objetivo de este ejercicio es obtener los límites departamentales a partir del shape de municipios de Colombia, este descargado de la página oficial del SIGOT. Este archivo puede ser de utilidad al momento de hacer el mapa de macrolocalización de algún departamento o municipio, en donde se muestre en este mapa, los departamentos y con algún otro color la zona dónde se encuentra el mapa.

Nota: el tutorial presente fue realizado con el software ArcGIS 10.3, de igual manera la herramienta se encuentra en versiones anteriores.

Primero abrimos el shapefile en ArcGIS usando el icono resaltado en amarillo:

Deberemos tener en el espacio de trabajo de ArcGIS algo como esto:

La herramienta de la cual se hará uso es “Dissolve”, esta se encuentra dentro del “ArcToolbox  Data Management Tools  Generalization  Dissolve”.

Se desplegara una ventana, en el archivo de entrada va a ir el municipio, en el campo “output” va el archivo de salida (es decir, donde queremos el nuevo shape que se creara), luego escogemos el campo sobre el cual se ejecutara el comando, para este caso escogemos el nombre del departamento, lo que significa que los polígonos (municipios) que hay dentro de cada departamento desaparecerán. Si quisiéramos solamente el límite de Colombia no escogemos ningún campo.

El archivo generado debe ser algo similar a lo siguiente:

El archivo usado en el presente tutorial se encuentra disponible en este link.

Extracción por mascara en R

Extracción por máscara en R

El objetivo del presente ejercicio es realizar la extracción por mascara (corte) de un grupo de archivos raster a partir de un shapefile, este usado como máscara, mediante un ciclo en R, es un proceso común cuando se tienen, por ejemplo, la precipitación (o alguna otra variable) para los 12 meses del año de todo Colombia y se quiere obtener solamente la precipitación de una unidad administrativa de menor área, para el presente caso es el departamento del Valle del Cauca; evitando hacer el proceso repetitivo, por ejemplo, en ArcGIS usando la herramienta “Extract by Mask”. Se hace uso de los ciclos, básicamente consiste en decirle al sistema que repita en cada una de las capas de un listado el corte a partir de una capa que representa el límite de nuestra área de estudio.
 Fuente: esri
Primero habilitamos las librerías de las cuales haremos uso:
require(raster) #si la librería no la tiene incluida se puede descargar usando install.packages(“raster”)
require(rgdal)
Luego abrimos el shapefile, que representa la máscara, bajo el comando “shapefile”.
mascara = shapefile(“E:/Blogger/Post_7/_data/_mask/valle_cacuca.shp”)
Para poder visualizar ambos archivos hacemos uso del comando “plot”.
plot(mascara)
class(mascara) #permite ver el tipo de archivo que es mascara, para este caso es un shapefile tipo polígono.
Luego hacemos el uso del comando “paste” para indicarle a R que lea todos los archivos que contienen la dirección escrita, esto permite unir la dirección, el 1 a 12 significa cada mes (es como finaliza el archivo) y por último la extensión que es .asc.
precipitacion=paste(“E:/Blogger/Post_7/_data/_prec_col/prec_”,1:12,”.asc”, sep=””)
Seguido de ello hacemos uso del comando “lapply” para indicar que todos los archivos son raster.
rasters=lapply(precipitacion, FUN = raster)
Luego realizamos el ciclo, aspectos a tener en cuenta:
·         El for es el comando que se utiliza
·         La “i” es para indicar que cada vez que esta letra cambie, conforme a la secuencia (1:12) será un archivo distinto al cual se le hace el proceso.
·         Lo que está dentro de las “{}” será lo que repetirá una y otra vez hasta finalizar.
·         El comando “crop” se usa para realizar el corte (equivalente al “extract by mask” de ArcGIS), en la primera posición va el dato de entrada, y en la segunda posición va el dato de máscara crop(rasters, mask),
·         El comando “mask” se usa para indicar que la máscara (extensión del archivo) será equivalente al objeto “mascara”.
·         El comando writeRaster es para guardar los archivos en una dirección de nuestro computador, y el comando “names” se utiliza para que guarde el archivo bajo el mismo nombre que el archivo de entrada.
for(i in 1:19){
  corte = crop(rasters[[i]],mascara)
  corte = mask(corte, mascara)
  writeRaster(corte,paste(“E:/Blogger/Post_7/_data/_prec_valle/”,names(rasters[[i]]),”.tif”,sep=””))
}
plot(corte[[1]])



De este link se puede descargar los datos que se hacen uso en este tutorial.
Gracias a @AntonioPantojaC por su colaboración!