#install.packages("knitr")
#require("knitr")#url = "http://geonetwork.humboldt.org.co/geonetwork/srv/api/records/e29b399c-24ee-4c16-b19c-be2eb1ce0aae/attachments/LHFI.zip"
#destfile = "D:/LHFI.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 IHEH para los años 1970, 1990, 2000 y 2015. Originalmente, el índice tiene un valor de 0 a 100 indicando en orden ascendente el grado de impacto humano. De acuerdo con el flujo de trabajo definido por el equipo de BioTablero, el IHEH se reclasifica para todos los años en cuatro categorías: 0-15 (Natural), 16-40 (Baja), 41-60 (Media) y 61-100 (Alta). Esta reclasificación se basa en Ayram et al. (2020). En primer lugar, la rutina lee cada uno de los archivos raster desde su ruta de ubicación con la función open de la librería rasterio de esta manera: rasterio.open(raster_fn), donde raster fn de corresponde a cada uno de los de los años del índice ya señalados 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 Huella.
- 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 IHEH, 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 IHEH 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:
# >>>IHEH = np.array([[12, 34, 68, 100], [25, 13, 48, 40],
# [40, 50, 40, 60], [0, 12, 3, 20]])
# >>> IHEH
# array([[12, 34, 68, 100],
# [25, 13, 48, 40],
# [40, 50, 40, 60],
# [0, 12, 3, 20]])
# >>> value_map = [((0, 15), 0), ((16, 40), 1),((41, 60), 2),((61, 100), 3)]
# >>> IHEH_cat=reclassify(arr, value_map)
# >>> IHEH_cat
# array([[0, 1, 3, 3],
# [1, 0, 2, 1],
# [1, 2, 1, 2],
# [0, 1, 0, 1]])Una vez se encuentra reclasificado el IHEH, 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. 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 IHEH reclasificado y convertirlo a un geodataframe.
- Operación: Toma el IHEH 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 IHEH 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:
# >>> IHEH_cat
# array([[0, 1, 3, 3],
# [1, 0, 2, 1],
# [1, 2, 1, 2],
# [0, 1, 0, 1]]), dtype=uint8)
# >>> IHEH_cat_shapes = rasterio.features.shapes(IHEH_cat)
# >>> IHEH_cat_shapes
# <generator object shapes at 0x7f9ace463820>
#
# >>> IHEH_cat_shape=shapes_to_geodataframe(IHEH_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 objeto geodataframe con los valores reclasificados del IHEH, se intersecta con la capa de geocercas Utilizando la función overlay de la librería geopandas, de acuerdo con el flujo de trabajo definido por el equipo de BioTablero. La capa de geocercas reúne todas las áreas de consulta por las que el usuario puede visualizar la información del IHEH y se debe llamar con anterioridad con el comando geofences = geopandas.read_file(geofences_path). La función opera de la siguiente manera:
- 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 geocercas acá llamada geofences y las opera en una intersección.
- Argumentos: Los dos primeros argumentos serán las capas a superponer de tipo shape (en este caso IHEH_cat_shape y geofences) y el tercero, 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: En el siguiente ejemplo se muestra la función y su operación:
# IHEH_cat_inters = intersection = geopandas.overlay(geofences, IHEH_cat_shape, how="intersection")
#intersection
# geofences IHEH_cat_shape 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 recategorizado y las geocercas en el anterior paso, se le agregan tres campos de acuerdo con el flujo de trabajo definido por el equipo de BioTablero: categoría (hf_cat), año (hf_year) y promedio (hf_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
# IHEH_cat_inters[hf_field_names.get("category")] = categories
#En donde category_map= {0: "natural", 1: "baja", 2: "media", 3: "alta"}
# IHEH_cat_inters[hf_field_names.get("category")] = categories
# Para el caso del campo del año:
# year = re.findall(r"(?:19|20)\d{2}", raster_fn)[0] #Se extre el año del nombre de la capa raster original del IHEH
# intersection[hf_field_names.get("year")] = year
#En donde hf_field_names = {
# "area": "area_ha",
# "average": "hf_avg",
# "category": "hf_cat",
# "persistence": "hf_pers",
# "protection": "binary_protected",
# "year": "hf_year"El campo del promedio del valor del índice para cada año por cada polígono dentro del geodataframe IHEH_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 IHEH a cada una de los años, polígonos y geocercas.
- Argumentos: Hace uso de una capa geodataframe producto de la intersección anterior (IHEH_cat_inters) y el archivo raster raster original del IHEH.
- Resultado: Genera el cálculo del promedio del IHEH a cada una de los años, polígonos y geocercas.
- Ejemplo:
#means = [item.get("mean") for item in zonal_stats(intersection, raster_fn)]
# intersection[hf_field_names.get("average")] = meansEn el paso 4.en el que se intersectan el IHEH (IHEH_cat_shape) y las geocercas (geofences), 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 generar uno único 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 IHEH_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 agregar 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 IHEH para cada geocerca intersectada.
- Resultado: Genera un nuevo geodataframe disuelto por los argumentos especificados.
- Ejemplo:
# Create empty GeoDataFrame to store the result.
# IHEH_final = geopandas.GeoDataFrame()
# IHEH_final = IHEH_final.append(intersection, ignore_index=True)
#dissolve_fields = list(geofences.columns) + hf_dissolve_fields
# dissolve_fields.remove("geometry")
# if None in dissolve_fields:
# dissolve_fields.remove(None)
# IHEH_final = IHEH_final.dissolve(by=dissolve_fields, aggfunc="mean")
# IHEH_final = IHEH_final.reset_index()
#Donde hf_dissolve_fields = [
# HF_FIELD_NAMES.get("category"),
# HF_FIELD_NAMES.get("year"),
# HF_FIELD_NAMES.get("protection")
#]Finalmente, la rutina calcula el área en hectáreas para cada uno de los polígonos que de la capa IHEH_final utilizando la función IHEH_final.geometry.get para obtener el área que se guarda implícitamente en algún lado del objeto geoDataFrame IHEH_final y conla función IHEH_final.geometry.arealo agrega al campo area anteriormente agregado; finalmente, el área es multiplicado por un facotr definido como 0.0001.
# Calcula el área solo si la capa está proyectada
# if IHEH_final.crs.is_projected:
# IHEH_final[hf_field_names.get("area")] = IHEH_final.geometry.area * area_factor
# else:
# warnings.warn("Area was not computed because input data is not projected.")
# Donde
#area_factor = 0.0001Referencias: Ayram, Camilo, Andrés Etter, Julián Díaz, Susana Rodríguez, Wilson Ramírez, and Germán Crozo. 2020. “Spatiotemporal Evaluation of the Human Footprint in Colombia : Four Decades of Anthropic Impact in Highly Biodiverse Ecosystems.” Ecological Indicators 117(April):106630.