<jupyter_start><jupyter_text>Indices E y CTakahashi, K., <NAME>., <NAME>., and <NAME> ( 2011), ENSO regimes: Reinterpreting the canonical and Modoki El Niño, Geophys. Res. Lett., 38, L10704, doi:10.1029/2011GL047364. El paper desarrollado por Takahashi et. al. propone que los dos primeros modos resultantes del análisis de funciones ortogonales empíricas (EOF) para las anomalías de temperatura superficial del mar no representan fenómenos distintos, EOF1 asociado a El Niño la Oscilación Sur (ENOS) y EOF2 asociado a El Niño Modoki, sino a la evolución no lineal del ENOS. Esta evolución se traduce a los índices E y C los cuales no guardan correlación (no se encuentran relacionados por construcción) y logran explicar la variabilidad presente en el Pacífico Ecuatorial, con el índice C representando el régimen de variabilidad en el Pacífico Central y el índice E en el Pacífico Este.El trabajo original hace uso del set de datos de temperatura superficial del mar de Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST) version 1.1; el trabajo de monitoreo por parte del IGP cuenta con el cálculo de este índice usando ERSSTv3b al cual se puede acceder mediante el siguiente [link](http://www.met.igp.gob.pe/variabclim/indices.html)Para realizar los cálculos usaremos los datos de temperatura ERSSTv5 [1] a los cuales accederemos mediante el DODS de IRI Data Library. Seguiremos los siguientes pasos:    Cargamos los datos usando xarray arreglando los tiempos de ser necesario restringiendo nuestro dominio a 10°S-10°N, 110°E-90°W y el tiempo desde enero de 1870 en adelanteCalculamos la anomalía tomando como periodo base 1979-2009. Este periodo será usado para todos los cálculos (EOFs y correlaciones)Realizamos el cálculo de EOFs sobre el periodo de la climatología.Proyectamos la serie completa de las anomalías sobre las EOFs calculadas en la climatología para obtener los componentes principales (PC). Esto se obtiene mediante un producto interno.Normalizamos los componentes principales usando la desviación estandar del periodo base y se aplica un filtro 1-2-1. Esto último lo logramos haciendo una correlación cruzadaRealizamos una regresión lineal de los campos de anomalía con los componentes principales para obtener los patrones espaciales de temperatura asociado a cada índice Referencias:[1]: <NAME>., <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, and <NAME>, 2017: Extended Reconstructed Sea Surface Temperature, version 5 (ERSSTv5): Upgrades, validations, and intercomparisons. J. Climate, 30, 8179-8205, doi:10.1175/JCLI-D-16-0836.1.  Implementación en PythonHaremos uso principalmente de dos librerias: xarray para el manejo de datos y eofs para realizar el calculo de EOFs sobre objetos de xarray. Los gráficos los manejaremos con matplotlib y cartopy.<jupyter_code>import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

import xarray as xr
import numpy as np
from eofs.xarray import Eof
from scipy.stats import linregress

%matplotlib inline<jupyter_output><empty_output><jupyter_text>1. Cargamos los datos directo de IRI Library<jupyter_code>ersstv5_iri = xr.open_dataset('http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCDC/.ERSST/.version5/.sst/dods',decode_times=False)
ersstv5_iri.T.attrs['calendar'] = '360_day'

# Ahora le decimos a xarray que interprete los tiempos
ersstv5_iri = xr.decode_cf(ersstv5_iri).rename({'X':'lon','Y':'lat','T':'time'}).isel(zlev=0).drop('zlev')
# Convertimos los tiempos a objetos de pandas
ersstv5_iri['time'] = ersstv5_iri.indexes['time'].to_datetimeindex()
ersstv5_iri
from scipy import signal
def xdetrend(xrobj, dim):
    mask = np.isnan(xrobj)
    out = xr.apply_ufunc(signal.detrend, xrobj.interpolate_na(dim='time').fillna(0), input_core_dims=[[dim]], output_core_dims=[[dim]], kwargs={'axis':-1})
    out = out.where(~mask)
    return out.transpose('time','lat','lon')
sst_pac = ersstv5_iri.sst.sel(lat=slice(-30, 30), lon=slice(110, 290), time=slice('1870-01-01', None)).load()
sst_pac = xdetrend(sst_pac, 'time')
sst_pac<jupyter_output><empty_output><jupyter_text>2. Calculamos la climatología y la anomalía<jupyter_code>clim = sst_pac.sel(time=slice('1979-01-01','2009-12-30')).groupby('time.month').mean(dim='time')
sst_anom = sst_pac.groupby('time.month') - clim
sst_anom<jupyter_output><empty_output><jupyter_text>3. Calculamos las EOFs para el periodo de la climatología<jupyter_code>clim_period = sst_anom.sel(time=slice('1979-01-01','2009-12-31'), lat=slice(-10, 10))
coslat = np.cos(np.deg2rad(clim_period.lat.data))
wgts = np.sqrt(coslat)[..., np.newaxis]

# Creamos un objeto solver que realizara los calculos de EOFs
# y le agregamos pesos a los datos de acuerdo a la raiz cuadrada
# del coseno de la latitud
solver = Eof(clim_period, weights=wgts)
clim_eof = solver.eofsAsCovariance(neofs=2)
clim_var = solver.varianceFraction(3)
clim_eof
import hvplot.xarray
solver.reconstructedField(2).hvplot(groupby='time',x='lon',y='lat', cmap=plt.get_cmap('RdBu_r')).redim.range(sst=(-1.5,1.5))
crs = ccrs.PlateCarree()
fig, axs = plt.subplots(dpi=200,
                        figsize=(8,3),
                        nrows=2,
                        subplot_kw={'projection':ccrs.PlateCarree(central_longitude=180)})
clim_eof.isel(mode=0).plot(ax=axs[0], transform=crs, add_colorbar=False)
clim_eof.isel(mode=1).plot(ax=axs[1], transform=crs, add_colorbar=False)
for ax in axs:
    ax.add_feature(cfeature.COASTLINE)<jupyter_output><empty_output><jupyter_text>Vemos que el modo 2 se encuentra invertido por lo que creamos un factor de corrección<jupyter_code>corr_factor = xr.DataArray(np.array([1,-1]),coords=[('mode',[0,1])])
clim_eof = clim_eof * corr_factor
clim_eof
crs = ccrs.PlateCarree()
fig, axs = plt.subplots(dpi=200,
                        figsize=(8,3),
                        nrows=2,
                        subplot_kw={'projection':ccrs.PlateCarree(central_longitude=180)})
clim_eof.isel(mode=0).plot(ax=axs[0], transform=crs, add_colorbar=False)
clim_eof.isel(mode=1).plot(ax=axs[1], transform=crs, add_colorbar=False)
for ax in axs:
    ax.add_feature(cfeature.COASTLINE)<jupyter_output><empty_output><jupyter_text>4. Proyectamos las anomalías sobre las EOFs de nuestro periodo base<jupyter_code>proj_sst = solver.projectField(sst_anom.sel(lat=slice(-10, 10)), neofs=2) * corr_factor
proj_sst<jupyter_output><empty_output><jupyter_text>5. Normalizamos y aplicamos filtro 1-2-1<jupyter_code>def ax_correlate(*args, **kwargs):
    return np.apply_along_axis(np.correlate, -1, *args, **kwargs)

def xsmooth(xrobj, weights):
    return xr.apply_ufunc(ax_correlate, xrobj, weights, input_core_dims=[['time'],[]], output_core_dims=[['time']], kwargs={'mode':'same'})
climpc_std = solver.pcs(npcs=2).std(dim='time')

# Normalizamos por la desviación estandar del periodo base
proj_sst_norm = proj_sst/climpc_std
PC1 = proj_sst_norm.sel(mode=0)
PC2 = proj_sst_norm.sel(mode=1)
C = (PC1+PC2)/(2**0.5)
E = (PC1-PC2)/(2**0.5)
fig, ax = plt.subplots()
E.sel(time=slice('1990-01-01',None)).plot(ax=ax)
ax.set_ylim(-2,5)
fig, ax = plt.subplots()
C.sel(time=slice('1990-01-01',None)).plot(ax=ax)
ax.set_ylim(-4,3)<jupyter_output><empty_output><jupyter_text>Podemos comparar nuestros resultados con los obtenidos por el IGP en su web de [Indices climáticos](http://www.met.igp.gob.pe/variabclim/indices.html)![eyc](http://www.met.igp.gob.pe/datos/EC.gif)<jupyter_code>style = dict(lw=0.5,c='k')
E_plot = E.sel(time=slice('1988-01-01', None))
C_plot = C.sel(time=slice('1988-01-01', None))

fig, ax = plt.subplots(nrows=2, dpi=200)

E_plot.plot(ax=ax[0], **style)
ax[0].set_ylim([-2, 5])

C_plot.plot(ax=ax[1], **style)
ax[1].set_ylim([-4, 3])

fig.subplots_adjust(hspace=0.6)<jupyter_output><empty_output><jupyter_text>5. Regresion lineal anom - pcs<jupyter_code>import tqdm.notebook as tqdm

def get_lincoeff(pc, anom, **lbar_kwargs):
    lin_coeff = np.zeros_like(anom[0])
    for i in tqdm.trange(anom.shape[1],**lbar_kwargs):
        for j in range(anom.shape[2]):
            lin_coeff[i,j] = linregress(pc.data,anom[:,i,j].data)[0]
    lin_coeff = xr.DataArray(lin_coeff, dims=['latitude','longitude'], coords=[anom.lat.data,anom.lon.data])
    return lin_coeff
E_pattern = get_lincoeff(E.sel(time=clim_period.time), sst_anom.sel(time=clim_period.time), desc='E pattern')
C_pattern = get_lincoeff(C.sel(time=clim_period.time), sst_anom.sel(time=clim_period.time), desc='C pattern')
C_pattern.plot.contourf(levels=np.arange(-1.3,1.31,0.2))
E_pattern.plot.contourf(levels=np.arange(-1.3,1.31,0.2))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12725789535019227, 'fraction_numerical': 0.034494814124227946, 'mean_word_length': 5.61171032357473, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '192211', 'license': 'MIT', 'n_tokens_mistral': 3488, 'n_tokens_neox': 3313, 'n_words': 761, 'path': 'Aplicaciones/IndiceEC.ipynb', 'repo_name': 'DangoMelon/python-geociencias', 'size': '8647'}