Curso: Procesamiento de datos geoespaciales de biodiversidad mediante el lenguaje de programación Python

Proyecto Final.

Alberto Méndez Rodríguez

Script realizado tomando como base el material desarrollado por Manuel Vargas, profesor curso: Procesamiento de datos geoespaciales de biodiversidad mediante el lenguaje de programación Python con la colaboración y auspicio de la ** Red de Ciencia de Datos para la Conservación de la Biodiversidad Mesoamericana ** (Red Bioma)

El enlace del script base es: [Material Curso] (https://datos-geoespaciales-biodiversidad.github.io/python/#procesamiento-de-datos-geoespaciales-de-biodiversidad-mediante-el-lenguaje-de-programacion-python)

Se escogieron las especies Cedro (Cedrela odorata), Almendro (Dipteryx oleifera) y Cristobal (Platymiscium parviflorum) por se especies en algunos casos protegidas pero todas importantes para la construcción de desarrollo de souveniers desde hace mucho tiempo en el Costa Rica.

Madera de las especies seleccionadas (en su orden): Cedrela odorata, Dipteryx oleifera, Platymiscium parviflorum

La información de avistamientos (presencia) se dividio por provincia para ver su distribución dentro del territorio nacional

# @title Instalacion algunas librerias
# Instalación de pygbif
!pip install pygbif --quiet
# Biblioteca requerida para mapas interactivos
!pip install mapclassify --quiet
# Biblioteca rasterio
!pip install rasterio --quiet
! pip install tabulate --quiet
!pip install folium --quiet

#@title Importar librerias
from pygbif import occurrences
import pandas as pd # biblioteca para análisis de datos
import matplotlib.pyplot as plt
import geopandas as gpd
import seaborn as sns
import plotly.express as px
import mapclassify
import folium
import rasterio
from folium.plugins import MarkerCluster
import folium
import rasterio
from folium import raster_layers
import tabulate

#@title Direccion archivo del google drive
from google.colab import drive
drive.mount('/content/drive')
Provincia_dire="/content/drive/My Drive/Colab Notebooks/datos/Provincias_WGS.shp"

Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount("/content/drive", force_remount=True).

#@title Lista de especies (solo 3)
lista_especie=["Cedrela odorata","Dipteryx oleifera","Platymiscium parviflorum"]

#@title Busqueda especie 0
# Nombre científico de la especie a buscar
especie = lista_especie[0]
limite = 300   # Límite de registros por solicitud
offset = 0    # Desplazamiento para iniciar la solicitud
registros_acumulados = []  # Lista para acumular los resultados
# Ciclo que obtiene los registros en la cantidad especificada
# en la variable limite
while True:
    # Solicitud
    res = occurrences.search(
        scientificName=especie,
        country	= 'CR',
        hasCoordinate=True,
        hasGeospatialIssue=False,
        limit=limite,
        offset=offset
    )

    # Extraer resultados
    registros = res.get("results", [])

    # Si ya no hay resultados, se detiene el ciclo
    if not registros:
        break

    # Agregar registros nuevos al acumulado
    registros_acumulados.extend(registros)

    # Actualizar el offset para la siguiente "página"
    offset += limite

# Convertir todo a un DataFrame de pandas
presencia0 = pd.DataFrame(registros_acumulados)

# Cantidad de registros de presencia recuperados
print(f"Se recuperaron {len(presencia0)} registros de presencia de {especie}")

Se recuperaron 143 registros de presencia de Cedrela odorata

#@title Busqueda especie 1
# Nombre científico de la especie a buscar
especie = lista_especie[1]
limite = 300   # Límite de registros por solicitud
offset = 0    # Desplazamiento para iniciar la solicitud
registros_acumulados = []  # Lista para acumular los resultados
# Ciclo que obtiene los registros en la cantidad especificada
# en la variable limite
while True:
    # Solicitud
    res = occurrences.search(
        scientificName=especie,
        country	= 'CR',
        hasCoordinate=True,
        hasGeospatialIssue=False,
        limit=limite,
        offset=offset
    )

    # Extraer resultados
    registros = res.get("results", [])

    # Si ya no hay resultados, se detiene el ciclo
    if not registros:
        break

    # Agregar registros nuevos al acumulado
    registros_acumulados.extend(registros)

    # Actualizar el offset para la siguiente "página"
    offset += limite

# Convertir todo a un DataFrame de pandas
presencia1 = pd.DataFrame(registros_acumulados)

# Cantidad de registros de presencia recuperados
print(f"Se recuperaron {len(presencia1)} registros de presencia de {especie}")

Se recuperaron 95 registros de presencia de Dipteryx oleifera

#@title Busqueda especie 2
# Nombre científico de la especie a buscar
especie = lista_especie[2]
limite = 300   # Límite de registros por solicitud
offset = 0    # Desplazamiento para iniciar la solicitud
registros_acumulados = []  # Lista para acumular los resultados
# Ciclo que obtiene los registros en la cantidad especificada
# en la variable limite
while True:
    # Solicitud
    res = occurrences.search(
        scientificName=especie,
        country	= 'CR',
        hasCoordinate=True,
        hasGeospatialIssue=False,
        limit=limite,
        offset=offset
    )

    # Extraer resultados
    registros = res.get("results", [])

    # Si ya no hay resultados, se detiene el ciclo
    if not registros:
        break

    # Agregar registros nuevos al acumulado
    registros_acumulados.extend(registros)

    # Actualizar el offset para la siguiente "página"
    offset += limite

# Convertir todo a un DataFrame de pandas
presencia2 = pd.DataFrame(registros_acumulados)

# Cantidad de registros de presencia recuperados
print(f"Se recuperaron {len(presencia2)} registros de presencia de {especie}")

Se recuperaron 76 registros de presencia de Platymiscium parviflorum

#@title Union de presencia de los 3 archivos

presencia_acumulada=pd.concat([presencia0,presencia1,presencia2], axis=0)
print(f"Se recuperaron {len(presencia_acumulada)} registros de presencia")
Presencia=presencia_acumulada[['species', 'decimalLongitude', 'decimalLatitude']]

Se recuperaron 314 registros de presencia

#@title Crear un GeoDataFrame a partir del DataFrame
Presencia_gdf = gpd.GeoDataFrame(
    Presencia,
    geometry=gpd.points_from_xy(presencia_acumulada.decimalLongitude, presencia_acumulada.decimalLatitude),
    crs='EPSG:4326'
)

#@title Cargar provincias de CR
Provincias_gdf = gpd.read_file(Provincia_dire)
Provincias_gdf

	PROVINCIA	PROV	geometry
0	ALAJUELA	2.0	POLYGON ((-84.66645 11.07237, -84.60605 11.038...
1	CARTAGO	3.0	POLYGON ((-83.81696 10.08363, -83.77386 10.062...
2	GUANACASTE	5.0	MULTIPOLYGON (((-85.22978 9.73273, -85.23029 9...
...	...	...	...
4	LIMON	7.0	MULTIPOLYGON (((-82.61599 9.63098, -82.6157 9....
5	PUNTARENAS	6.0	MULTIPOLYGON (((-82.96208 8.21685, -82.96227 8...
6	SAN JOSE	1.0	POLYGON ((-83.93499 10.10305, -83.93492 10.101...

7 rows × 3 columns

#@title Unión (join) espacial
especies_provincias_gdf = Presencia_gdf.sjoin(
    Provincias_gdf,
    predicate='intersects'
)

#@title Conteo de avistamientos de especies por provincia
conteo_especies= especies_provincias_gdf['PROVINCIA'].value_counts() # especies_provincias_gdf.groupby("PROVINCIA").species.nunique()
# Renombrar todas las columnas
#conteo_especies.columns = ['Provincia', 'Provincia']
conteo_especies = conteo_especies.rename_axis('Provincia').reset_index(name='Count')
conteo_especies.columns = ['PROVINCIA', 'Cantidad']
conteo_especies

print("Tabla de Avistamientos de las 3 especies seleccionadas")
display(conteo_especies.style)

Tabla de Avistamientos de las 3 especies seleccionadas

	PROVINCIA	Cantidad
0	GUANACASTE	95
1	HEREDIA	73
2	SAN JOSE	45
3	PUNTARENAS	39
4	LIMON	34
5	ALAJUELA	21
6	CARTAGO	3

#@title Union bases de Provincias y conteo especies
provincia_union=pd.merge(Provincias_gdf,conteo_especies,on='PROVINCIA',how='inner')

# @title Gráfico de avistamientos por especie
# 1. Contar registros por país
registros_x_especie = especies_provincias_gdf['species'].value_counts()
bar_labels=registros_x_especie.index
bar_colors = ['purple', 'blue', 'green']

# 2. Generar un gráfico de barras
plt.figure(figsize=(10, 6)) # tamaño del gráfico
plt.bar(registros_x_especie.index, registros_x_especie.values,label=bar_labels,color=bar_colors)

plt.legend(loc='upper right',title='Especies',)
plt.xticks(rotation=60)

# 3. Especificar el título y las etiquetas del gráfico
plt.title("Figura 1. \nCantidad de Avistamientos por especie \npara las 3 especies seleccionadas\n",
          fontdict={'fontsize': 18, 'color':"blue",'fontweight': 'bold'})
plt.xlabel('Especie')
plt.ylabel('Cantidad de registros')

# 4. Mostrar el gráfico
plt.show()

Se observa que la especie Cedrela odorata 3s la que tiene mayor cantidad de avistamientos, la razón de esta es que no es una especie en peligro de extinción, aunque si es madera de muy buena calidad

#@title Grafico interactivo de barras apiladas de avistamientos por especie y provincia
species_provincia_counts = especies_provincias_gdf.groupby(['species', 'PROVINCIA'])['species'].count().reset_index(name='count')

# 2. Crear el gráfico de barras con Plotly Express (cambiando x e y)
fig = px.bar(species_provincia_counts,
             x='PROVINCIA',  # Provincia en el eje x
             y='count',
             color='species',  # Especie como color
             title='Figura 2. Cantidad de avistamientos por provincia y especie',
             labels={'PROVINCIA': 'Provincia', 'count': 'Cantidad de registros', 'species': 'Especie'})

# 3. Mostrar el gráfico
fig.show()

Se puede observar en la figura anterior que el Cedro tiene una distribución muy parecida en Guanacaste, Heredia, Puntarenas y San José, pues es una especie muy adaptable. Mientras el Cristobal ha sido mas observado en Guanacaste, este comportamiento se puede deber principalmente que aunque existe en otras provincias por ejemplo Puntarenas se encuentra más en bosques densos, mientras en Guanacaste se pueden observar como arboles aislados

#@title Grafico interactivo Distribucion de avistamientos por provincia y especie
species_provincia_counts = especies_provincias_gdf.groupby(['PROVINCIA', 'species'])['species'].count().reset_index(name='count')

# Crear el histograma
fig = px.histogram(species_provincia_counts,
                   x="species",
                   y="count",
                   color="PROVINCIA",
                   text_auto=True,
                   title="Figura 3. Distribución de avistamientos por especies por provincia",
                   labels={'PROVINCIA': 'Provincia', 'count': 'Cantidad de registros', 'species': 'Especie'},
                   barmode='group') # Apilar las barras para cada especie

fig.show()

Se puede ratificar lo observado en la Figura 2, pues el Cristobal tiene una gran cantidad de avistamientos en Guanacaste.El Cedro la distribución es parecida en las provinvias de Guanacaste, Heredia, Puntarenas y San José.El caso del Almendro los avistamientos mayores son el la provincia de Heredia, donde se ubica la Estación La Selva, que es una zona muy visitada y el almendro además es muy frecuente en esa zona en especial hacia el norte de la provincia.

#@title Mapa Folium de avistamientos divididos por especie a nivel nacional
mapa=[]

popup = folium.GeoJsonPopup(
    fields=["species"],
    localize=True,
    labels=True,
    style="background-color: yellow;",
)

tooltip = folium.GeoJsonTooltip(
    fields=["species"],
    localize=True,
    sticky=False,
    labels=True,
    style="""
        background-color: #F0EFEF;
        border: 2px solid black;
        border-radius: 3px;
        box-shadow: 3px;
    """,
    max_width=800,
)
title_html = '''
             <div style="position: absolute; top: 10px; left: 50%; transform: translate(-50%, 0%); z-index: 9999; font-size:20px; background-color: white; padding: 5px;">
                 <strong>Mapa de avistamientos por especie</strong>
             </div>
             '''
mapa=folium.Map(location=[9.8,-84.2],zoom_start=8)
mapa.get_root().html.add_child(folium.Element(title_html))

#folium.Map(location=[9.8,-84.2],zoom_start=10).add_to(mapa)

folium.GeoJson(Provincias_gdf,
               name="Provincias",
               color="red",
               weight=0.7,
               fill=False,
               show=True,
               ).add_to(mapa)

# definir color
color_dict = {
    "Cedrela odorata": "purple",
    "Dipteryx oleifera": "blue",
    "Platymiscium parviflorum": "green"
}

#puntos por especie
for index, row in especies_provincias_gdf.iterrows():
    folium.CircleMarker(
        location=[row['decimalLatitude'], row['decimalLongitude']],
        radius=3,
        color=color_dict.get(row['species'], 'red'),  # Get color from dictionary, default to red
        fill=True,
        fill_color=color_dict.get(row['species'], 'red'),
        fill_opacity=0.7,
         tooltip=row['species']  # Display species name in popup
    ).add_to(mapa)


folium.LayerControl().add_to(mapa)


leyenda_html = """
     <div style="position: fixed;
                 bottom: 50px; left: 50px; width: 200px; height: 90px;
                 border:2px solid grey; z-index:9999; font-size:14px;
                 background-color:white;
                 ">
         <div style="text-align: center; font-weight: bold; padding-top: 5px;">Leyenda</div>
         <p style="margin: 2px;"><i class="fa fa-circle" style="color:purple"></i> Cedrela odorata</p>
         <p style="margin: 2px;"><i class="fa fa-circle" style="color:blue"></i> Dipteryx oleifera</p>
         <p style="margin: 2px;"><i class="fa fa-circle" style="color:green"></i> Platymiscium parviflorum</p>
     </div>
     """
mapa.get_root().html.add_child(folium.Element(leyenda_html))
mapa

Make this Notebook Trusted to load map: File -> Trust Notebook

#@title Mapa Folium de avistamientos agrupados de las 3 especies.
mapa=[]

popup = folium.GeoJsonPopup(
    fields=["species"],
    localize=True,
    labels=True,
    style="background-color: yellow;",
)

tooltip = folium.GeoJsonTooltip(
    fields=["species"],
    localize=True,
    sticky=False,
    labels=True,
    style="""
        background-color: #F0EFEF;
        border: 2px solid black;
        border-radius: 3px;
        box-shadow: 3px;
    """,
    max_width=800,
)
title_html = '''
             <div style="position: absolute; top: 10px; left: 50%; transform: translate(-50%, 0%); z-index: 9999; font-size:20px; background-color: white; padding: 5px;">
                 <strong>Mapa de presencia de especies agrupadas</strong>
             </div>
             '''
mapa=folium.Map(location=[9.8,-84.2],zoom_start=8)
mapa.get_root().html.add_child(folium.Element(title_html))

#folium.Map(location=[9.8,-84.2],zoom_start=10).add_to(mapa)

folium.GeoJson(Provincias_gdf,
               name="Provincias",
               color="red",
               weight=0.7,
               fill=False,
               show=True,
               ).add_to(mapa)
marker_cluster = MarkerCluster(name="Reportes presencia de las especies").add_to(mapa)
folium.GeoJson(especies_provincias_gdf,
               popup=popup,
               tooltip=tooltip,
               marker=folium.CircleMarker(radius=5, color="red", fill_color="red")
               ).add_to(marker_cluster)

folium.LayerControl().add_to(mapa)

mapa

Make this Notebook Trusted to load map: File -> Trust Notebook

#@title Coropleta
mapa_coropletas=[]
title_html = '''
             <div style="position: absolute; bottom: 10px; left: 30%; transform: translate(-50%, 0%); z-index: 9999; font-size:20px; background-color: white; padding: 5px;">
                 <strong>Mapa de Coropletas de avistamientos </br> de las 3 especies por Provincia</strong>
             </div>
             '''
mapa_coropletas=folium.Map(location=[9.8,-84.2],zoom_start=8)

mapa_coropletas.get_root().html.add_child(folium.Element(title_html))

folium.Choropleth(
    geo_data=provincia_union,
    name="choropleth",
    data=provincia_union,
    columns=["PROVINCIA", "Cantidad"],
    key_on="feature.properties.PROVINCIA",  # Make sure this matches the actual name in your GeoJSON
    fill_color="viridis",
    fill_opacity=0.7,
    line_opacity=0.2,
    legend_name="Cantidad de Especies"
).add_to(mapa_coropletas)

mapa_coropletas

Make this Notebook Trusted to load map: File -> Trust Notebook

#@title Llamar imagen Sentinel
imagen_sentinel = rasterio.open(
    'https://github.com/datos-geoespaciales-biodiversidad/python/raw/refs/heads/main/datos/sentinel/gandoca-20240105.tif'
)

#@title Ver imagen Sentinal en color verdadero
red_band = imagen_sentinel.read(4).astype('float32')  # Assuming band 4 is red
green_band = imagen_sentinel.read(3).astype('float32') # Assuming band 3 is green
blue_band = imagen_sentinel.read(2).astype('float32') # Assuming band 2 is blue


# Valores nulos reemplazados con minimo valor de la banda
red_band = np.nan_to_num(red_band, nan=red_band.min())
green_band = np.nan_to_num(green_band, nan=green_band.min())
blue_band = np.nan_to_num(blue_band, nan=blue_band.min())

# Normalizar bandas de 0-1
red_norm = (red_band - red_band.min()) / (red_band.max() - red_band.min())
green_norm = (green_band - green_band.min()) / (green_band.max() - green_band.min())
blue_norm = (blue_band - blue_band.min()) / (blue_band.max() - blue_band.min())

# eliminar valores fuera rango normalizacion
red_norm = red_norm.clip(0, 1)
green_norm = green_norm.clip(0, 1)
blue_norm = blue_norm.clip(0, 1)

# Agrupar bandas normalizadas en RGB
rgb = np.dstack((red_norm, green_norm, blue_norm))

# Ver la imagen RGB
plt.figure(figsize=(8, 8))
plt.imshow(rgb)
plt.title("Imagen Sentinel - RGB")
plt.axis('off')
plt.show()

#@title Imagen Sentinel aclarada
red_band = imagen_sentinel.read(4).astype('float32')  # se asume banda 4 como roja
green_band = imagen_sentinel.read(3).astype('float32') # se asume banda 3 como verde
blue_band = imagen_sentinel.read(2).astype('float32') # se asume banda 2 como azul


# reemplazo valores nulos con minimo valor de banda
red_band = np.nan_to_num(red_band, nan=red_band.min())
green_band = np.nan_to_num(green_band, nan=green_band.min())
blue_band = np.nan_to_num(blue_band, nan=blue_band.min())

# Normalizar valores de bandas de 0 - 1
red_norm = (red_band - red_band.min()) / (red_band.max() - red_band.min())
green_norm = (green_band - green_band.min()) / (green_band.max() - green_band.min())
blue_norm = (blue_band - blue_band.min()) / (blue_band.max() - blue_band.min())

# eliminar los valores de banda fuera del rango
red_norm = red_norm.clip(0, 1)
green_norm = green_norm.clip(0, 1)
blue_norm = blue_norm.clip(0, 1)

# se incrementa brillo para mas claridad
brightness_factor = 1.5  # se incrementa brillo un 50%
red_norm = np.clip(red_norm * brightness_factor, 0, 1)
green_norm = np.clip(green_norm * brightness_factor, 0, 1)
blue_norm = np.clip(blue_norm * brightness_factor, 0, 1)


# Agrupar bandas ara RGB
rgb = np.dstack((red_norm, green_norm, blue_norm))

# Ver imagen RGB clara
plt.figure(figsize=(8, 8))
plt.imshow(rgb)
plt.title("Imagen Sentinel - RGB (Brillo aumentado 50%)")
plt.axis('off')
plt.show()

#@title Grabar imagen
plt.savefig('rgb_image.jpg')

<Figure size 1250x1000 with 0 Axes>