Monitor de sequía de México situación de las cuencas del río Yaqui y río Sonora

Por Luis Alan Navarro Navarro, 24 mayo 2022, Lnavarro@colson.edu.mx

¿Por qué estas dos cuencas?

Primero, mi foco de interés es el abastecimiento de agua de la ciudad de Hermosillo. Las fuentes de agua, superficiales o subterráneas, provienen de escurrimientos o inflitraciones a los acuíferos dentro de estas cuencas hidrológicas. La ciudad en suma, depende de la lluvia que se precipita en estas dos cuencas.

Fuentes de agua o captaciones de la ciudad

Son de dos tipos: Agua de pozo o agua tomada de una presa. Las presas son tres: Abelardo L. Rodríguez, El Molinito y El Novillo; de estas el agua se extrae para ser potabilizada. Sin embargo, requieren cierto nivel de almacenamiento para que se pueda extraer el agua óptimamente: Abelardo L. Rodríguez, 40% de su capacidad; El Molinito, 60% de su capacidad; y El Novillo al 30% de su capacidad. Segundo, la capacidad de las potabilizadoras tiene un límite: 3,300 lps (para transformar el flujo a litros anuales: lps x 60 s. x 60 min. x 24 hr. x 365 días o los días que opera).

Por último, la recarga de los pozos también proviene del agua de lluvia que se infiltra en el subsuelo. La ciencia geológica ha delimitado estos reservorios subterráneos de agua y los llama acuíferos. La mayoría de los pozos están en el acuífero Mesa del Serí-La Victoria (al oriente de la ciudad); luego en el acuífero Costa de Hermosillo (al poniente de la ciudad); también están los acuíferos del Río San Miguel (al norte) y La Poza (al sur).

"Bottom line"

Dependemos de la lluvia. Las presas y los reservorios subterráneos (acuíferos) solo dan amortiguamiento a los años de lluvias escasas.

¿Qué sería no depender de las lluvias? Por ejemplo, desalinizar el agua de mar; o reusar las aguas negras tratadas (que en cierta forma dependen de la fuente de agua potable).

El código

Como primer paso descarga dos archivos en ESRI-shp creados por el autor en base a la división municipal de INEGI y las cuencas delineadas por la CONAGUA. Guarda los archivos en tu carpeta o directorio de trabajo.

Descargar aquí.

Monitor de sequía

Una buena herramienta para monitorear la evolución de la sequía es el Monitor de sequía de México del Sistema Meteorológico Nacional, CONAGUA. Visita esta página. Y descarga archivo de municipios (botón azul en medio de la página). Es un archivo de una hoja de cálculo en formato MS Excel o "csv".

Disclaimer

No soy experto en Python, por lo que entiendo que este código es perfectible. Es más, invito a los usuarios a que perfeccionen estos mapas y el código.

### módulos
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
import geopandas as gpd
import contextily as cx
import pandas as pd
import os
import re
import fiona
import numpy as np
import contextily as ctx

# Paso 1. El archivo del monitor de sequía
# define la ruta a tu directorio de trabajo
Ruta0_main = '/Users/COLSON/Desktop/Anaconda_files/Indicadores_HCV_Agua'
os.chdir(Ruta0_main)
## Excel
sequiaAll = pd.read_excel('MunicipiosSequia-3.xlsx',index_col=None, sheet_name='MUNICIPIOS')

# Paso 2. Cambio el directorio de trabajo 
Ruta_rioSon = '/Users/COLSON/Desktop/Anaconda_files/Indicadores_HCV_Agua/municipios_RS_RY'
os.chdir(Ruta_rioSon)

# Paso 3. Cargo el primero shapefile, el del río Sonora
fp = "municipios_RS_UTM.shp"
    ## Read file using gpd.read_file()
mpios_rson = gpd.read_file(fp)

   ## Extraemos el listado de municipios
mpios_son = mpios_rson.loc[:,'CVEGEO']
mpios_son = mpios_son.tolist()

# Paso 4. Ahora el río Yaqui
fp = "municipios_RY_UTM.shp"
    ## SHAPEFILE DEL RIO YAQUI
    ## Read file using gpd.read_file()
mpios_yaqui = gpd.read_file(fp)

    ## Extraemos el listado de municipios del río Yaqui
mpios_son_yaqui = mpios_yaqui.loc[:,'CVEGEO']
mpios_son_yaqui = mpios_son_yaqui.tolist()

# Paso 5 Uno ambas listas y evito duplicados
list_1 = mpios_son_yaqui
set_1 = set(mpios_son_yaqui)
set_2 = set(mpios_son)

list_2_items_not_in_list_1 = list(set_2 - set_1)
combined_list = list_1 + list_2_items_not_in_list_1
    ## Hay que filtrar la base de datos Sonora y Chihuahua primero
sequia_son = sequiaAll[sequiaAll['CVE_CONCATENADA'].isin(combined_list)]

# Paso 6. Pasamos a la forma larga
sequia_sonL = sequia_son.melt(id_vars= ['CVE_CONCATENADA', 'CVE_ENT','CVE_MUN','NOMBRE_MUN', \
                          'ENTIDAD','ORG_CUENCA','CLV_OC','CON_CUENCA','CVE_CONC'],
                       value_vars=sequia_son.columns[9:sequia_son.shape[1]], var_name = 'quincena' , 
                        value_name ='sequia' )

# Paso 7. AQUI SELECCIONAMOS LA FECHA QUE NOS INTERESA
    # Mayo 15 del 2022, define un rango en el que cae esta fecha 
mask = (sequia_sonL['quincena'] > '2022-05-10') & (sequia_sonL['quincena'] <= '2022-05-16')
    ################################################       
    # print(sequia_sonL.loc[mask])
quincena_sequia = sequia_sonL.loc[mask]

#### verifico que no esté vacía la selección
quincena_sequia 

	CVE_CONCATENADA	CVE_ENT	CVE_MUN	NOMBRE_MUN	ENTIDAD	ORG_CUENCA	CLV_OC	CON_CUENCA	CVE_CONC	quincena	sequia
23170	8009	8	9	Bocoyna	Chihuahua	Río Bravo	VI	Rio Bravo	12	2022-05-15	D2
23171	8012	8	12	Carichí	Chihuahua	Río Bravo	VI	Rio Bravo	12	2022-05-15	D2
23172	8013	8	13	Casas Grandes	Chihuahua	Río Bravo	VI	Rio Bravo	12	2022-05-15	D2
23173	8017	8	17	Cuauhtémoc	Chihuahua	Río Bravo	VI	Rio Bravo	12	2022-05-15	D2
23174	8018	8	18	Cusihuiriachi	Chihuahua	Río Bravo	VI	Rio Bravo	12	2022-05-15	D2
...	...	...	...	...	...	...	...	...	...	...	...
23235	26067	26	67	Villa Hidalgo	Sonora	Noroeste	II	Rios Yaqui y Matape	4	2022-05-15	D2
23236	26068	26	68	Villa Pesqueira	Sonora	Noroeste	II	Rios Yaqui y Matape	4	2022-05-15	D2
23237	26069	26	69	Yécora	Sonora	Noroeste	II	Rios Yaqui y Matape	4	2022-05-15	D2
23238	26071	26	71	Benito Juárez	Sonora	Noroeste	II	Rios Yaqui y Matape	4	2022-05-15	D2
23239	26072	26	72	San Ignacio Río Muerto	Sonora	Noroeste	II	Rios Yaqui y Matape	4	2022-05-15	D2

70 rows × 11 columns

# Paso 7
## Create a dictionary where you assign each attribute value to a particular color
color_mapping = {"D0": "#FDF733", "D1": "#FBD37F", "D2": "#E5982E", "D3": "#E8442E","D4": "#741D12", 
                 "No_S.": "#d6e250" }

# Paso 8
## DE AQUI SE DERIVA EL ANÁLISIS PARA UNA U OTRA CUENCA para esta quincena específica
quincena_sequia = quincena_sequia.rename(columns={"CVE_CONCATENADA": "CVEGEO"})

# Paso 9 río Sonora
#######################################
    ## como le agrego a cada shp la variable de sequía???
    ## https://stackoverflow.com/a/40603298
mpios_rson['seq1'] = mpios_rson.CVEGEO.map(quincena_sequia.set_index('CVEGEO')['sequia'].to_dict())

quincena_sequia['CVEGEO']= quincena_sequia['CVEGEO'].astype(str)
mpios_rson['seq1'] = mpios_rson['CVEGEO'].map(quincena_sequia.set_index('CVEGEO')['sequia'])
mpios_rson.loc[mpios_rson['seq1'].isnull(), 'seq1'] = "No_S."

# Paso 10. comenzamos con la grafica para el rio Sonora
sequia_rson_por = mpios_rson.groupby('seq1')['area_por'].agg(['sum'])

# Paso 11 checamos que categorias tenemos 
sequia_rson_por

	sum
seq1
D2	0.8217
D3	0.1783

# Categorias:
names  = {"D0": 'Anormalmente seco',
           "D1": 'Sequía moderada',
           "D2": 'Sequía severa',
           "D3": 'Sequía extrema',
           "D4": 'Sequía excepcional',
           "No_S.": 'No sequía',
          }

# Paso 12 esto se hace manualmente en fucnión del Paso 11
  
    ## solo para las dos clases presentes
# D0_v = round(sequia_rson_por.at['D0','sum']*100,2)
# D1_v = round(sequia_rson_por.at['D1','sum']*100,2)
D2_v = round(sequia_rson_por.at['D2','sum']*100,2)
D3_v = round(sequia_rson_por.at['D3','sum']*100,2)
# No_S_v = round(sequia_rson_por.at['No_S.','sum']*100,2)

# Paso 12 HACEMOS LA LEYENDA DEL MAPA manualmente

    ## Patch(facecolor="#FDF733", edgecolor='black',
    ##                         label= 'D0 Anormalmente seco '+ str(D0_v)+'% de la cuenca.'),
    ## Patch(facecolor="#d6e250", edgecolor='black',
    ##                         label='No_S. No sequía '+ str(No_S_v)+'% de la cuenca.')
    ## Patch(facecolor="#FBD37F", edgecolor='black',
    ##                         label='D1 Sequía moderada '+ str(D1_v1)+'% de la cuenca.'),
    ## Patch(facecolor="#E8442E", edgecolor='black',
    ##                         label='D3 Sequía extrema '+ str(D3_v1)+'% de la cuenca.')

    ##                Patch(facecolor="#FDF733", edgecolor='black',
    ##                     label= 'D0 Anormalmente seco '+ str(D0_v2)+'% \n de la cuenca.'),
    ##                Patch(facecolor="#FBD37F", edgecolor='black',
    ##                     label='D1 Sequía moderada '+ str(D1_v2)+'%  \n de la cuenca.'),                    
    ##                Patch(facecolor="#E5982E", edgecolor='black',
    ##                     label='D2 Sequía severa '+ str(D2_v2)+'% de la cuenca.'),
    ##                Patch(facecolor="#d6e250", edgecolor='black',
    ##                    label='No_S. No sequía '+ str(No_S_v2)+'% de la cuenca.')                 
                  
legend_elements2 = [ 
                       Patch(facecolor="#E5982E", edgecolor='black',
                            label='D2 Sequía severa '+ str(D2_v)+'% de la cuenca.'),
                        Patch(facecolor="#E8442E", edgecolor='black',
                            label='D3 Sequía extrema '+ str(D3_v)+'% de la cuenca.')
                  ]

# Paso 13 Buscamos otro background para el mapa
    ## to the Pseudo-Mercator Projection
mpios_rson_PSM = mpios_rson.to_crs(epsg=3857)

# Paso 14 ajustamos la geometría al nuevo sistema de coordenadas
mpios_rson_PSM['coords'] = mpios_rson_PSM['geometry'].apply(lambda x: x.representative_point().coords[:])
mpios_rson_PSM['coords'] = [coords[0] for coords in mpios_rson_PSM['coords']]

El mapa se crea también en alta resolución (300 DPI) y se descarga a tu directorio de trabajo.

# Control figure size in here
fig, ax = plt.subplots(figsize=(15,15))
# Plot the data
# Loop through each attribute type and plot it using the colors assigned in the dictionary
for ctype, data in mpios_rson_PSM.groupby('seq1'):    
    # Define the color for each group using the dictionary
    color = color_mapping[ctype]
    label = names[ctype]
    # Plot each group using the color defined above
    data.plot(color=color,
              ax=ax,
              edgecolor="black",
              legend=True,
              alpha=0.5)
ax.set_title('Monitor de Sequía de México, cuenca río Sonora, mayo 15 2022',fontweight="bold", size=14) # Title
# ax.set(title='Monitor de Sequía de México, cuenca río Sonora, octubre 15 2020')
for idx, row in mpios_rson_PSM.iterrows():
   plt.annotate(s=row['seq1'], xy=row['coords'], horizontalalignment='center', color='black')
ax.legend(handles=legend_elements2, loc="upper left", fontsize=12)
# Add basemap
ctx.add_basemap(ax)
fig.savefig('fig_mayo15_2022.png', dpi = 300)

<ipython-input-22-3ee56544299f>:18: MatplotlibDeprecationWarning: The 's' parameter of annotate() has been renamed 'text' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
  plt.annotate(s=row['seq1'], xy=row['coords'], horizontalalignment='center', color='black')

Para la cuenca del río Yaqui

    ## 1 código compacto
mpios_yaqui['seq1'] = mpios_yaqui['CVEGEO'].map(quincena_sequia.set_index('CVEGEO')['sequia'])
mpios_yaqui.loc[mpios_yaqui['seq1'].isnull(), 'seq1'] = "No_S."
y2022 = mpios_yaqui.groupby('seq1')['area_por'].agg(['sum'])
print(y2022)

          sum
seq1         
D1     0.0174
D2     0.6208
D3     0.0994
No_S.  0.2619

    ## 2 Creo etiquetas
# D0_v2 = round(y101521.at['D0','sum']*100,2)
D1_v2 = round(y2022.at['D1','sum']*100,2)
D2_v2 = round(y2022.at['D2','sum']*100,2)
D3_v2 = round(y2022.at['D3','sum']*100,2)
No_S_v2 = round(y2022.at['No_S.','sum']*100,2)
    
    ## 3 Creamos centroides para etiquetas
mpios_yaqui['coords'] = mpios_yaqui['geometry'].apply(lambda x: x.representative_point().coords[:])
mpios_yaqui['coords'] = [coords[0] for coords in mpios_yaqui['coords']]
    ## 4 AQUI CAMBIAMOS DE PROYECCION PARA USAR DIFERENTE BACKGROUND
    ## to the Pseudo-Mercator Projection
mpios_yaqui_PSM = mpios_yaqui.to_crs(epsg=3857)
    ##
    ## 5 VUELVO A CREAR LOS CENTROIDES
mpios_yaqui_PSM['coords'] = mpios_yaqui_PSM['geometry'].apply(lambda x: x.representative_point().coords[:])
mpios_yaqui_PSM['coords'] = [coords[0] for coords in mpios_yaqui_PSM['coords']]
    ##
    ## 6 HACEMOS LA LEYENDA DEL MAPA

    ## Patch(facecolor="#FDF733", edgecolor='black',
    ##                         label= 'D0 Anormalmente seco '+ str(D0_v)+'% de la cuenca.'),
    ## Patch(facecolor="#d6e250", edgecolor='black',
    ##                         label='No_S. No sequía '+ str(No_S_v)+'% de la cuenca.')
    ## Patch(facecolor="#FBD37F", edgecolor='black',
    ##                         label='D1 Sequía moderada '+ str(D1_v1)+'% de la cuenca.'),
    ## Patch(facecolor="#E8442E", edgecolor='black', 
    #                 label='D3 Sequía extrema '+ str(D3_v1)+'% de la cuenca.')

legend_elements3 = [ 
                    Patch(facecolor="#FBD37F", edgecolor='black',
                         label='D1 Sequía moderada '+ str(D1_v2)+'%  \n de la cuenca.'),                    
                    Patch(facecolor="#E5982E", edgecolor='black',
                         label='D2 Sequía severa '+ str(D2_v2)+'% de la cuenca.'),
                    Patch(facecolor="#E8442E", edgecolor='black', 
                         label='D3 Sequía extrema '+ str(D3_v2)+'% de la cuenca.'),
                    Patch(facecolor="#d6e250", edgecolor='black',
                        label='No_S. No sequía '+ str(No_S_v2)+'% de la cuenca.')
                  ]

    ## 7 
    ## Control figure size in here
fig, ax = plt.subplots(figsize=(15,15))

    ## Plot the data

    ## Loop through each attribute type and plot it using the colors assigned in the dictionary
for ctype, data in mpios_yaqui_PSM.groupby('seq1'):    
    ## Define the color for each group using the dictionary
        color = color_mapping[ctype]
        label = names[ctype]
    ## Plot each group using the color defined above
        data.plot(color=color,
              ax=ax,
              edgecolor="black",
              legend=True,
              alpha=0.5)

    ax.set_title('Monitor de Sequía de México, cuenca río Yaqui, mayo 15 2022',fontweight="bold", size=14) # Title
    ## ax.set(title='Monitor de Sequía de México, cuenca río Sonora, octubre 15 2021')
for idx, row in mpios_yaqui_PSM.iterrows():
   plt.annotate(s=row['seq1'], xy=row['coords'], horizontalalignment='center', color='black')
ax.legend(handles=legend_elements3, loc="upper right", fontsize=12)
    ## Add basemap
ctx.add_basemap(ax, zoom=9)
fig.savefig('figYaqui22.png', dpi = 200)

<ipython-input-26-250b918382fb>:22: MatplotlibDeprecationWarning: The 's' parameter of annotate() has been renamed 'text' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
  plt.annotate(s=row['seq1'], xy=row['coords'], horizontalalignment='center', color='black')

Muchas gracias por tu atención.