proyecciones_hist.py

# -*- coding: utf-8 -*-
"""
Editor de Spyder

Este es un archivo temporal.
"""

####################################################################################################
##### Importamos las bibliotecas necesarias para la adquisición de datos y manejo de variables #####
####################################################################################################
import netCDF4 as nc4
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate
import cartopy.feature as cfeat
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import shapefile
from shapely.geometry import shape, Point

historico = ['/home/meteo/datos/CMIP5/CCSM4/sud_tas_Amon_CCSM4_1pctCO2_r1i1p1_194801-200512.nc',
             '/home/meteo/datos/CMIP6/EC_Earth3/sesa_tas_Amon_EC-Earth3_historical_r1i1p1f1_gr_194801-201412.nc',
             '/home/meteo/datos/CMIP6/EC_Earth3-Veg/sesa_tas_Amon_EC-Earth3-Veg_historical_r1i1p1f1_gr_194801-201412.nc',
             '/home/meteo/datos/CMIP6/CAMS-CSM1-0/sesa_tas_Amon_CAMS-CSM1-0_historical_r1i1p1f1_gn_194801-201412.nc',
             '/home/meteo/datos/CMIP6/CanESM5_CCCma/sud_tas_Amon_CanESM5_historical_r1i1p1f1_gn_195001-201412.nc',
             '/home/meteo/datos/CMIP6/BCC_CSM2/sud_tas_Amon_BCC-CSM2-MR_historical_r1i1p1f1_gn_195001-201412.nc',
             '/home/meteo/datos/CMIP6/MIROC6/sud_tas_Amon_MIROC6_historical_r1i1p1f1_gn_195001-201412.nc',
             '/home/meteo/datos/CMIP6/MRI_ESM2/sud_tas_Amon_MRI-ESM2-0_historical_r1i1p1f1_gn_195001-201412.nc']
observado = '/home/meteo/datos/Uy/CRU_tmp_sud_1948_2017.nc'
modelos = ['CCSM4','EC_Earth3','EC_Earth3-Veg','CAMS','CanESM5','BCC_CSM2','MIROC6','MRI_ESM2']

####################################################################################################
############      Acomodamos los datos para que tengan el mismo rango temporal     #################
####################################################################################################
tas2 = []
lon = []
lat = []
for i in range(len(modelos)):
    ncin = nc4.Dataset(historico[i], 'r')
    tas2.append(ncin.variables['tas'][:])
    lon.append(ncin.variables['lon'][:])
    lat.append(ncin.variables['lat'][:])
    ncin.close()

tas = []
alcance = min(tas2[0].shape[0],tas2[1].shape[0],tas2[2].shape[0],tas2[3].shape[0])
for i in range(len(modelos)):
    if i < 4:
        tas.append(tas2[i][24:int(alcance),:,:]) # los datos de estos modelos están desde 1948
    else:
        tas.append(tas2[i][0:int(alcance)-24,:,:]) # los datos de estos modelos están desde 1950
        

#cargamos las variables del observado tmp
ncin = nc4.Dataset(observado, 'r')
tmp = ncin.variables['tmp'][:]
lono = ncin.variables['lon'][:]
lato = ncin.variables['lat'][:]
ncin.close()
tmp=tmp[24:int(alcance),:,:]
tmp[tmp > 1000] = np.nan

####################################################################################################
###############     Definimos funciones para utilizar en el resto del código      ##################
####################################################################################################
def regrillado(obs,mod,lonobs,latobs,lonmod,latmod):
    '''función para regrillar las observaciones (resolución alta) a la grilla del modelo (resolución más baja) para una comparación apropiada'''
    #old grid dim
    loni=lonobs           # 'obs' son las temperaturas observadas CRU
    lati=latobs
    #new grid dim
    lonn=lonmod           # 'mod' son las temperaturas del modelo CMIP
    latn=latmod
    #create mesh
    X, Y = np.meshgrid(loni, lati)
    XI, YI = np.meshgrid(lonn,latn)
    X = X + 360          # esto es porque las observaciones tienen valores de longitud en (-180,180]
    #interp
    A2 = obs.transpose([1,2,0]).reshape(int(len(loni)*len(lati)),mod.shape[0])
    new_A=scipy.interpolate.griddata((X.flatten(),Y.flatten()),A2,(XI,YI),method='linear')
    new_A = new_A.transpose([2,0,1])
    return new_A    # devuelve el arreglo original regrillado según lonmod y latmod.

def climatologia(temp,latitud,longitud):
    '''función para calcular la climatología de temp(tiempo,lon,lat)'''
    new_temp = temp.reshape(int(temp.shape[0]/12),12,int(latitud.shape[0]),int(longitud.shape[0]))
    return np.mean(new_temp, axis = 0)

def graf_clim(clim,lati,long,model,mes):
    '''función para graficar las climatologías de cada mes'''
    plt.figure()
    meses = ['enero','febrero','marzo','abril','mayo','junio','julio','agosto','setiembre','octubre','noviembre','diciembre']
    clevs = [0,3,6,9,12,15,18,21,24,27,30,33]
    #clevs = [-30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30]
    ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=-56))
    fill = ax.contourf(np.real(long), np.real(lati), np.real(clim), np.real(clevs), transform=ccrs.PlateCarree(), cmap=plt.cm.coolwarm, extend='both')
    ax.coastlines()
    #ax.gridlines()
    ax.add_feature(cfeat.RIVERS)
    ax.add_feature(cfeat.BORDERS)
    ax.set_extent([min(long), max(long), max(lati), min(lati)], crs=ccrs.PlateCarree())
    ax.set_xticks([295, 300, 305, 310], crs=ccrs.PlateCarree())
    ax.set_yticks([-25, -30, -35], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=True,
                                       number_format='.0f')
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    plt.title(model+' climatología '+ meses[mes])
    plt.colorbar(fill, orientation='horizontal')
    #plt.title('prec mm/día', fontsize=12)
    plt.savefig('/home/meteo/investigacion/WR.slp/proyecciones/'+model+'/temp/clim_'+model+'%d.png' % (mes), dpi=300)

def puntos_adentro_Uy(x,y,polygon):
    '''función para evaluar si una coordenada espacial está dentro del territorio uruguayo'''
#    r = shapefile.Reader('/home/meteo/investigacion/WR.slp/URY_adm/URY_adm0.shp')
    # get the shapes
#    shapes = r.shapes()
    # build a shapely polygon from your shape
#    polygon = shape(shapes[0])    
    # build a shapely point from your geopoint
    point = Point(x, y)
        # the contains function does exactly what you want
    return polygon.contains(point)

def dif_clim_tri(obs,mod,latitud,longitud,model,trimestre):
    '''función para calcular la diferencia de climatologías trimestrales'''
    # graficamos la diferencia de las climatologias trimestrales
    dif_clim = mod - obs
    plt.figure()
    clevs = [-6,-5,-4,-3,-2,-1, 0,1,2,3,4,5,6]
    ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=-56))
    fill = ax.contourf(np.real(longitud), np.real(latitud), np.real(dif_clim), np.real(clevs), transform=ccrs.PlateCarree(), cmap=plt.cm.coolwarm, extend='both')
    ax.coastlines()
    #ax.gridlines()
    ax.add_feature(cfeat.RIVERS)
    ax.add_feature(cfeat.BORDERS)
    ax.set_extent([min(longitud), max(longitud), max(latitud), min(latitud)], crs=ccrs.PlateCarree())
    ax.set_xticks([295, 300, 305, 310], crs=ccrs.PlateCarree())
    ax.set_yticks([-25, -30, -35], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=True,
                                       number_format='.0f')
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    plt.title(trimestre +' ('+model+' - CRU)')
    plt.colorbar(fill, orientation='horizontal')
    plt.savefig('/home/meteo/investigacion/WR.slp/proyecciones/'+model+'/temp/dif_clim_'+trimestre+'.png', dpi=300)
    
####################################################################################################
#######   Calculamos las climatologías mensuales de las observaciones y de los modelos   ###########
####################################################################################################
# calculamos la climatología mensual de los datos observados (CRU)
clim_cru = np.zeros((12,len(lato),len(lono)))
clim_cru = climatologia(tmp,lato,lono)
for m in range(12):
    graf_clim(clim_cru[m,:,:],lato,lono,'CRU',m)

# calculamos la climatología mensual para los modelos (CMIP5/6)
clim_modelos = []
for i in range(len(modelos)):
    clim_modelos.append(climatologia(tas[i] - 273.15,lat[i],lon[i]))
    for m in range(12):
        graf_clim(clim_modelos[i][m,:,:],lat[i],lon[i],modelos[i],m)

####################################################################################################
#######   Comparamos las climatologías trimestrales entre los modelos y las observaciones    #######
####################################################################################################
# Para una comparación más apropiada, regrillamos las observaciones (mejor resolución) a la grilla de cada modelo.
for i in range(len(modelos)):
    tmpr = regrillado(tmp,tas[i],lono,lato,lon[i],lat[i])   # REGRILLADO
    tmpr2 = tmpr.reshape(int(tmpr.shape[0]/12),12,int(lat[i].shape[0]),int(lon[i].shape[0]))

    tas3 = np.zeros((int(tas[i].shape[0]/12),12,int(lat[i].shape[0]),int(lon[i].shape[0])))
    tas3 = tas[i].reshape(int(tas[i].shape[0]/12),12,int(lat[i].shape[0]),int(lon[i].shape[0])) - 273.15

    # climatologia trimestrales de datos CRU regrillados
    tmpr_mam = np.mean(np.squeeze(np.mean(tmpr2[:,2:5,:,:], axis=1)),axis=0) 
    tmpr_jja = np.mean(np.squeeze(np.mean(tmpr2[:,5:8,:,:], axis=1)),axis=0) 
    tmpr_son = np.mean(np.squeeze(np.mean(tmpr2[:,8:11,:,:], axis=1)),axis=0) 
    tmpr_def = np.mean(np.squeeze(np.mean(tmpr2[:,[0,1,11],:,:], axis=1)),axis=0) 

    # climatologia trimestrales de datos CMIP5
    tas_mam = np.mean(np.squeeze(np.mean(tas3[:,2:5,:,:], axis=1)),axis=0)
    tas_jja = np.mean(np.squeeze(np.mean(tas3[:,5:8,:,:], axis=1)),axis=0)
    tas_son = np.mean(np.squeeze(np.mean(tas3[:,8:11,:,:], axis=1)),axis=0)
    tas_def = np.mean(np.squeeze(np.mean(tas3[:,[0,1,11],:,:], axis=1)),axis=0)

    # calculamos las diferencias y graficamos
    dif_clim_tri(tmpr_mam,tas_mam,lat[i],lon[i],modelos[i],'MAM')
    dif_clim_tri(tmpr_jja,tas_jja,lat[i],lon[i],modelos[i],'JJA')
    dif_clim_tri(tmpr_son,tas_son,lat[i],lon[i],modelos[i],'SON')
    dif_clim_tri(tmpr_def,tas_def,lat[i],lon[i],modelos[i],'DEF')

####################################################################################################
#####   Climatología de las obs y los modelos considerando un promedio espacial dentro de UY   #####
####################################################################################################
r = shapefile.Reader('/home/meteo/investigacion/WR.slp/URY_adm/URY_adm0.shp')   # carga el shapefile de Uruguay
shapes = r.shapes()                                                             # get the shapes
polygon = shape(shapes[0])                                                      # build a shapely polygon from your shape
clim_cru_uy = np.zeros((12))    # observaciones
for m in range(12):
    cont = 0
    for i in range(len(lato)):
        for j in range(len(lono)):
            if puntos_adentro_Uy(lono[j],lato[i],polygon):   # función definida arriba, evalúa si lono[j],lato[i] está dentro de Uy
                clim_cru_uy[m] += clim_cru[m,i,j]
                cont += 1
    print(m,cont)
    clim_cru_uy[m] /= cont
    
clim_mod_uy = np.zeros((8,12))    # modelos
for k in range(len(modelos)):
    for m in range(12):
        cont = 0
        for i in range(len(lat[k])):
            for j in range(len(lon[k])):
                if puntos_adentro_Uy(lon[k][j]-360,lat[k][i],polygon):
                    clim_mod_uy[k,m] += clim_modelos[k][m,i,j]
                    cont += 1
        clim_mod_uy[k,m] /= cont
    print(k,cont)
# graficamos las climatologías
plt.plot(clim_cru_uy, label = 'CRU', color= 'black')
plt.plot(clim_mod_uy[0,:], label = modelos[0])
plt.plot(clim_mod_uy[1,:], label = modelos[1])
plt.plot(clim_mod_uy[2,:], label = modelos[2])
plt.plot(clim_mod_uy[3,:], label = modelos[3])
plt.plot(clim_mod_uy[4,:], label = modelos[4])
plt.plot(clim_mod_uy[5,:], label = modelos[5])
plt.plot(clim_mod_uy[6,:], label = modelos[6])
plt.plot(clim_mod_uy[7,:], label = modelos[7])
plt.title('Climatología temperatura media Uruguay (1950 - 2005)')
plt.ylabel('$^\circ $C')
plt.legend(fancybox=True, framealpha=0.4)
plt.ylim((5,35))
plt.xticks(ticks=np.arange(12), labels=('Ene', 'Feb', 'Mar', 'Abr', 'May', 'Jun', 'Jul', 'Ago', 'Set', 'Oct', 'Nov', 'Dic'))
plt.savefig('/home/meteo/investigacion/WR.slp/proyecciones/clim_tas_hist.png', dpi=150)

####################################################################################################
########            Calculamos medidas para evaluar la habilidad de los modelos          ###########
####################################################################################################

# Calculamos la media de la climatología para las observaciones y para los modelos
media_uy_obs = np.mean(clim_cru_uy)
media_uy_mod = np.mean(clim_mod_uy, axis = 1)
sesgo_anual = np.zeros((8)) # Xm - Xobs anual
for i in range(8):
    sesgo_anual[i] = media_uy_mod[i] - media_uy_obs

# RMSE
rmse = np.zeros((8))
for i in range(8):
    for t in range(12):
        rmse[i] += (clim_mod_uy[i,t] - clim_cru_uy[t])**2
    rmse[i] = np.sqrt(rmse[i]/12)
rmse[rmse > 100] = np.nan

# calculo indice VI: es el cociente entre la desviación estándar de cada modelo y la desviación estándar de las observaciones
std_cru = np.std(tmp, axis = 0)
std_mod = []
for i in range(len(modelos)):
    std_mod.append(np.std(tas[i],axis = 0))

std_cru_uy = 0    # observaciones
cont = 0
for i in range(len(lato)):
    for j in range(len(lono)):
        if puntos_adentro_Uy(lono[j],lato[i],polygon):
            std_cru_uy += std_cru[i,j]
            cont += 1
std_cru_uy /= cont
    
std_mod_uy = np.zeros((8))    # modelos
for k in range(len(modelos)):
    cont = 0
    for i in range(len(lat[k])):
        for j in range(len(lon[k])):
            if puntos_adentro_Uy(lon[k][j] - 360,lat[k][i],polygon):
                std_mod_uy[k] += std_mod[k][i,j]
                cont += 1
    std_mod_uy[k] /= cont

VI = np.zeros((8))
VI = std_mod_uy/std_cru_uy