El siguiente documento tiene como fin describir y explicar las rutinas de código que permiten descargar el producto de Índice de Condición Estructural (SCI por su siglas en inglés) (Hansen et al. 2019) y posterior edición para su incorporación en el flujo de trabajo de la plataforma BioTablero. En primer lugar se detallan las funciones que descargan y descomprimen el producto a trabajar, posteriormente se toman como base las rutinas realizadas por Marcerlo Villa y Jaime Burbano relacionado con la edición de la capa de Índice espacial de la Huella Humana. El código base completo puede ser encontrado en la dirección de la plataforma Github: https://github.com/PEM-Humboldt/huella-humana-indicadores/tree/python. Igualmente, todos los argumentos creados para que las funciones se ejecuten de manera correcta se pueden encontrar en: https://github.com/PEM-Humboldt/huella-humana-indicadores/blob/python/src/utils/constants.py.
#install.packages("knitr")
#require("knitr")#url = "https://springernature.figshare.com/ndownloader/files/15667565"
#destfile = "D:/SCI.zip"#download.file(url = url, destfile = destfile) #Recuerde cambiar en este espacio con la ruta donde desea guardar el archivo descargado en su ordenador#unzip(zipfile = destfile, exdir = "D:/") #Recuerde verificar tanto la ruta del archivo, #como la dirección donde desea descomprimir la información. Dentro del archivo descargado y descomprimido en el paso anterior, se encuentra el Índice de Condición Estructural y trae clasificado de 1 a 18, donde 0 es un bajo y 18 un alto valor de integridad de los bosques. De acuerdo con el flujo de trabajo definido por el equipo de BioTablero, el SCI se reclasifica en dos categorías de la siguiente manera: 1-13 (Baja) y 14-18(Moderada).En primer lugar, la rutina lee el raster desde su ruta de ubicación con la función open de la librería rasterio de esta manera: rasterio.open(raster_sci), donde raster_sci corresponde al raster del índice ya descargado en formato .tif. Posteriormente, la reclasificación se hace con la función reclassify que tiene las siguientes características:
- Objetivo:Reclasificar el raster del Índice de Condición Estructural.
- Operación: La función toma los valores antiguos del raster y los reemplaza por los rangos de valores nuevos que el usuario requiera.
- Argumentos: Requiere de arr, un arreglo tipo numpy (matriz) en el que se almacenan los valores de cada uno de los archivos raster del SCI, y el argumento value_map, en el que se lista en un tipo tupla el límite inferior y superior de cada rango a reclasificar, y el nuevo valor que toma cada rango.
- Resultado: Genera un nuevo raster clasificado con los nuevos valores de categoría.
- Ejemplo: En el siguiente ejemplo se muestra el argumento SCI que indica el raster de entrada con los valores antiguos, y en el argumento valupe_map que contiene los rangos y nuevos valores del raster:
# >>>SCI = np.array([[12, 3, 8, 10], [15, 13, 18, 1],
# [1, 3, 4, 16], [0, 12, 3, 2]])
# >>> SCI
# array([[12, 3, 8, 10],
# [15, 13, 18, 1],
# [1, 3, 4, 16],
# [0, 12, 3, 2]])
# >>> value_map = [((0, 13), 0), ((14, 18), 1)]
# >>> SCI_cat=reclassify(arr, value_map)
# >>> SCI_cat
# array([[0, 0, 0, 0],
# [1, 0, 1, 0],
# [0, 0, 0, 1],
# [0, 0, 0, 0]])Una vez se encuentra reclasificado el SCI, se vectoriza en dos pasos. En primer lugar se convierte en un objeto espacial tipo shapes con la función shapes de la librería rasterio módulo features (rasterio.features.shapes). Posteriormente el reusltado se convierte en un objeto geodataframe con la función shapes_to_geodataframe que tiene las siguientes características:
- Objetivo: Convertir el raster de SCI reclasificado y convertilo a un geodataframe.
- Operación: Toma el SCI reclasificado en el paso anterior y convierte todas las características (o features) en un geoDataFrame.
- Argumentos: shapes que es de tipo shapefile y es el producto de la función rasterio.features.shapes, crs que será el sistema de referencia de la nueva capa a transformar en tipo texto, y field_name que es el nombre que llevará el campo dentro del objeto geodataframe y donde se almacenará el valor reclasificado del SCI en tipo texto.
- Resultado: Genera un nuevo objeto tipo geoDataFrame.
- Ejemplo: En el siguiente ejemplo se muestra en primer lugar la transformación del raster clasificado con la función rasterio.features.shapes y su posterior transformación a un geodataframe:
#>>> SCI_cat
# array([[0, 0, 0, 0],
# [1, 0, 1, 0],
# [0, 0, 0, 1],
# [0, 0, 0, 0]]), dtype=uint8)
# >>> SCI_cat_shapes = rasterio.features.shapes(SCI_cat)
# >>> SCI_cat_shapes
# <generator object shapes at 0x7f9ace463820>
#
# >>> SCI_cat_shape=shapes_to_geodataframe(SCI_cat_shapes, src.crs.to_string())
# value_field = features.columns[0]
# value geometry
# 0 1.0 POLYGON ((0.00000 0.00000, 0.00000 1.00000, 1....
# 1 1.0 POLYGON ((4.00000 0.00000, 4.00000 2.00000, 5....
# 2 0.0 POLYGON ((0.00000 1.00000, 0.00000 2.00000, 1....
# 3 1.0 POLYGON ((0.00000 2.00000, 0.00000 3.00000, 1....
# 4 1.0 POLYGON ((3.00000 2.00000, 3.00000 5.00000, 4....
# 5 0.0 POLYGON ((4.00000 2.00000, 4.00000 5.00000, 5....
# 6 0.0 POLYGON ((3.00000 0.00000, 3.00000 1.00000, 2....
# 7 1.0 POLYGON ((1.00000 4.00000, 1.00000 5.00000, 2....Una vez se haya creado el geodataframe del producto reclasificado del Índice de Condición Estructural, se intersecta con el producto de Huella Espacial Humana que se extrae de la rutina de Github de este mismo producto (la rutina del índice de huella se entrega junto con esta rutina con el nombre Rutina_IHEH_github.Rmd). El producto del Índice de Huella Humana además trae dentro de sus campos las geocercas y tipos de áreas protegidas por las cuales el usuario filtra la información dentro del Biotablero, de esta manera y con esta intersección, el Índice de Condición Estructural queda con estos campos agregados sin necesidad de ser intersectada con las geocercas que establece el instituto. La función que realiza esta intersección es overlay de la librería geopandas y tiene las siguientes características:
- Objetivo: Hacer una intersección entre dos capas tipo geodataframe.
- Operación: Toma el geoDataFrame generado en el anterior paso y la capa de la Huella Humana que es igualmente un geoDataFrame y las opera en una intersección.
- Argumentos: Los dos primeros argumentos serán las capas a superponer de tipo shape, en este caso SCI_final que es la capa final y procesada de Huella Humana y SCI_cat_shape. El tercer argumento es how, la forma en que las capas se quieren superponer (intersection, union, symetrical difference, difference) en tipo texto.
- Resultado: Genera un nuevo objeto tipo geodataframe.
- Ejemplo:
# SCI_cat_inters = intersection = geopandas.overlay(SCI_cat_shape, IHEH_final, how="intersection")
#intersection
# SCI_cat_shape IHEH_final geometry
#0 1 1 POLYGON ((1.000000000 2.000000000, 2.000000000...
#1 2 1 POLYGON ((2.000000000 2.000000000, 2.000000000...
#2 2 2 POLYGON ((3.000000000 4.000000000, 4.000000000...Al objeto geodataframe producto de la intersección entre las capas del IHEH y el SCI categorizado en el anterior paso, se le agregan tres campos de acuerdo con el flujo de trabajo definido por el equipo de BioTablero: categoría (sci_cat), año (sci_year) y promedio (sci_avg). Para lo anterior los nombres y abreviaciones se introducen en un vector que contiene los nombres de los campos de la siguiente manera:
# Para el caso de agregar el campo de categoría según la reclasificación:
# categories = intersection[value_field].map(category_map) #Se extraen los valores de la categoría de la capa intersectada
# SCI_cat_inters[sci_field_names.get("category")] = categories
#En donde category_map= {0: "baja", 1: "moderada"}
# SCI_cat_inters[sci_field_names.get("category")] = categories
# Para el caso del campo del año:
# year = re.findall(r"(?:19|20)\d{2}", raster_sci)[0] #Se extre el año del nombre de la capa raster original del SCI
# intersection[sci_field_names.get("period")] = year
#En donde sci_field_names = {
# "area": "area_ha",
# "average": "sci_avg",
# "category": "sci_cat",
# "protection": "binary_protected",
# "year": "sci_year"El campo del promedio del valor del índice para cada año por cada polígono dentro del geodataframe SCI_cat_inters y cada geocerca intersectada, se agrega con la función zonal_stats que tiene las siguientes características:
- Objetivo: Genera el cáculo de estadísticas zonales.
- Operación: Toma el geoDataFrame generado en el anterior paso de la intersección y calcula el valor promedio del SCI para cada geocerca.
- Argumentos: Hace uso de una capa geodataframe producto de la intersección anterior (SCI_cat_inters) y el archivo raster original del SCI.
- Resultado: Genera el cálculo del promedio del SCI a cada una de las geocercas.
- Ejemplo:
#means = [item.get("mean") for item in zonal_stats(intersection, raster_sci)]
# intersection[sci_field_names.get("average")] = meansEn el paso 4.en el que se intersectan la Huella Humana (IHEH_final) y el SCI categorizado y vectorizado (SCI_cat_inters), se generan 15 campos (“anu”, dcs, dnmi, drmi, pnn, pnr, rfpn, rfpr, rn, rnsc, sfa, sff, sfl, vp, ar) correspondientes a cada tipo de áreas protegidas existentes en el Registro Único Nacional de Áreas Protegidas (RUNAP). Para simplificar estos campos es necesario genear uno único campo que resume la presencia cada tipo de área protegida de acuerdo con el flujo de trabajo definido por el equipo de BioTablero, lo cual se hace uso de la función compute_protection_sequence. La función compute_protection_sequence tiene las siguientes características:
- Objetivo: Genera campo binario que da cuenta de la ausencia o presencia de área protegida.
- Operación: Asigna un valor de 0 (ausencia) o 1 (presencia) cuando el polígono correspondiente incluye o no área protegida. De acuerdo con el flujo definido por el equipo de BioTablero, cada uno de los campos de tipo de área protegida almacena el identificador único (ID) asignado para cada área protegida desde la capa original del RUNAP. La función convierte este campo a un valor de presencia (1) cuando el ID del área protegida es mayor a cero. Por defecto, la capa de geocercas codifica la ausencia de un área protegida en cada polígono con un valor de cero (0).
- Argumentos:df que se refiere a un dataframe (tabla) que almacena la secuencia de valores tomados de los 15 campos de áreas protegidas, y fields que es un vector que contiene el nombre de los campos a concatenar.
- Resultado: Genera un campo de tipo string, y corresponde a una cadena de 15 caracteres de acuerdo con la presencia o ausencia de cada uno de los tipos de áreas protegidas.
- Ejemplo:#>>> values = [[0, 0, 0, 27100002], [2060001, 0, 0, 0], [0, 0, 0,0]]
# >>> fields = ["anu","dcs","dnmi","drmi"]
# >>> df = pd.DataFrame(values, columns=fields)
# >>> df
# anu dcs dnmi drmi
# 0 0 0 0 27100002
# 1 2060001 0 0 0
# 2 0 0 0 0
# >>> compute_protection_sequence(df, fields=["anu","dnmi","drmi"])
# 0 001
# 1 100
# 2 000
# dtype: objectAl producto anterior se le aplica un dossolve a cada uno de los productos obtenidos con la función SCI_final.dissolve que tiene las siguientes características:
- Objetivo: Disolver geometrías.
- Operación: Toma el geoDataFrame generado en el anterior con los campos agregados y estadístiacas generadas.
- Argumentos: Hace uso del argumento by que se refiere a los campos por los cuales se quiere hacer el dissolve, en este caso será el año, la categoría del índice y el promedio del índice, en este caso un tipo lista. El segundo argumento es aggfunc que se refiere a la característica por la cual se requiere aregar el dissolve en este caso el promedio “mean” en tipo texto. En este caso se utiliza el promedio como función de agregación porque se quiere obtener el valor promedio del SCI para cada geocerca intersectada.
- Resultado: Genera un nuevo geodataframe disuelto por los argumentos especificados.
- Ejemplo:
# Create empty GeoDataFrame to store the result.
# SCI_final = geopandas.GeoDataFrame()
# SCI_final = SCI_final.append(intersection, ignore_index=True)
#dissolve_fields = list(geofences.columns) + sci_dissolve_fields
# dissolve_fields.remove("geometry")
# if None in dissolve_fields:
# dissolve_fields.remove(None)
# SCI_final = SCI_final.dissolve(by=dissolve_fields, aggfunc="mean")
# SCI_final = SCI_final.reset_index()
#Donde sci_dissolve_fields = [
# SCI_FIELD_NAMES.get("category"),
# SCI_FIELD_NAMES.get("year"),
# SCI_FIELD_NAMES.get("protection")
#]Finalmente, la rutina calcula el área en hecareas para cada uno de los polígonos que de la capa SCI_final utilizando la función SCI_final.geometry.get para obtener el área que se guarda implcitamente en algún lado del objeto geoDataFrame SCI_final y con la función SCI_final.geometry.arealo agrega al campo area anteriormente agregado; finalmente, el área es multiplicado por un factor definido como 0.0001.
# Calcula el área solo si la capa está proyectada
# if SCI_final.crs.is_projected:
# SCI_final[hf_field_names.get("area")] = SCI_final.geometry.area * area_factor
# else:
# warnings.warn("Area was not computed because input data is not projected.")
# Donde
#area_factor = 0.0001Referencias: Hansen, Andrew, Kevin Barnett, Patrick Jantz, Linda Phillips, Scott J. Goetz, Matt Hansen, Oscar Venter, James E. M. Watson, Patrick Burns, Scott Atkinson, Susana Rodríguez-Buritica, Jamison Ervin, Anne Virnig, Christina Supples, and Rafael De Camargo. 2019. “Global Humid Tropics Forest Structural Condition and Forest Structural Integrity Maps.” Scientific Data 6(1):232.