00:17:10	Teresa:	Buenas tardes con todo@s!
00:18:15	Christian Barreto:	buenas tardes
00:28:13	Christian Barreto:	import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import xarray as xr

dset = xr.open_dataset("precip.mon.mean.nc")
print(dset)
dset
00:29:17	Christian Barreto:	pp=dset.precip
data_c=pp.sel(time=slice("1981-01-01","2010-12-01"),lat=slice(10, -60),lon=slice(250, 350))

month_idxs=data_c.groupby('time.month').groups
# Extract the time indices corresponding to all the Januarys 
jan_idxs=month_idxs[1]

# Extract the january months by selecting 
# the relevant indices
data_jan=data_c.isel(time=jan_idxs)
data_jan
00:30:08	Christian Barreto:	#ONI
years_nino=[1982,1987,1991,1997,2002,2004,2009]
years_nina=[1985,1988,1998,1999,2000,2007,2010]
years_normal=[1981,1983,1984,1986,1989,1990,1992,1993,1994,1995,1996,2001,2003,2005,2006,2008]
00:30:44	Christian Barreto:	nino_idx=[]
nina_idx=[]
normal_idx=[]

year_idxs=data_jan.groupby("time.year").groups

for i in years_nino:
    nino_idx.append(year_idxs[i][0])
    
for i in years_nina:
    nina_idx.append(year_idxs[i][0])
    
for i in years_normal:
    normal_idx.append(year_idxs[i][0])
    
data_jan_nino=data_jan.isel(time=nino_idx)
data_jan_nina=data_jan.isel(time=nina_idx)
data_jan_normal=data_jan.isel(time=normal_idx)

data_jan_nino_mean=data_jan_nino.mean(dim="time",skipna=True)
data_jan_normal_mean=data_jan_normal.mean(dim="time",skipna=True)
00:31:05	Christian Barreto:	data_jan_nina_mean=data_jan_nina.mean(dim='time',skipna=True)
00:31:42	Christian Barreto:	import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib import colorbar, colors
import cartopy.feature as cf

# Draw coastlines of the Earth
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cf.BORDERS)
#ax.coastlines() 
ax.add_feature(cf.COASTLINE)
00:31:47	Christian Barreto:	#adding ejes
xticks=([-120,-100,-80,-60,-40,-20,0])
yticks=([-75,-60,-45,-30,-15,0,15])
ax.set_xticks(xticks, crs=ccrs.PlateCarree())
ax.set_yticks(yticks, crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,number_format='.1f')
lat_formatter = LatitudeFormatter(number_format='.1f')
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
00:31:53	Christian Barreto:	#adding grillas

ax.gridlines(draw_labels=False, xlocs=xticks, ylocs=yticks)
(data_jan_nino_mean-data_jan_normal_mean).plot(robust=True)

plt.title("nino-normal, 81-10")
ax.set_aspect('auto', adjustable=None)
00:35:22	Christian Barreto:	# con t-test
from scipy.stats import ttest_ind, ttest_ind_from_stats

data=data_jan_nino
pvalues = np.zeros((data.shape[1],data.shape[2]),dtype=float)

for ni in range(0,data.shape[2]):   # loop over longitudes
    for nj in range(0, data.shape[1]):  # loop over latitudes
00:35:28	Christian Barreto:	info_a=data_jan_nino.isel(lat=nj, lon=ni).values
        info_b=data_jan_normal.isel(lat=nj, lon=ni).values
        
        array_sum = np.sum(info_a)
        array_has_nan = np.isnan(array_sum)
          
        if array_has_nan == True:
00:35:34	Christian Barreto:	#            print("F")
            #r=np.nan
            #p=np.nan
        
            pvalues[nj,ni] = np.nan
            
        else:
#            print("ga")
                                 
            result = ttest_ind(info_a,info_b).pvalue

            pvalues[nj,ni] = result
00:36:09	Christian Barreto:	data_set = xr.Dataset( coords={'lon': ([ 'lon'], data.lon),
                                             'lat': (['lat',], data.lat)})
    
data_set["pvalues"] = (['lat', 'lon'],  pvalues)

data_set.to_netcdf("pvalues_nino_nnl_81-10.nc",mode='w')
00:36:27	Christian Barreto:	import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib import colorbar, colors
import cartopy.feature as cf

# Draw coastlines of the Earth
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cf.BORDERS)
#ax.coastlines() 
ax.add_feature(cf.COASTLINE)

#adding ejes
xticks=([-120,-100,-80,-60,-40,-20,0])
yticks=([-75,-60,-45,-30,-15,0,15])
ax.set_xticks(xticks, crs=ccrs.PlateCarree())
ax.set_yticks(yticks, crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,number_format='.1f')
lat_formatter = LatitudeFormatter(number_format='.1f')
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
00:36:35	Christian Barreto:	ax.gridlines(draw_labels=False, xlocs=xticks, ylocs=yticks)
00:36:38	Christian Barreto:	levels=[-2.5,-2,-1,0,1,2,3]
00:36:44	Christian Barreto:	data_set.pvalues.where(pvalues < 0.05).plot()
00:38:06	Christian Barreto:	plt.title('nino-normal 81-10')
plt.xlabel('longitud')
plt.ylabel('latitud')
00:38:09	Christian Barreto:	ax.set_aspect('auto', adjustable=None)
01:02:37	Christian Barreto:	https://docs.google.com/spreadsheets/d/1V2KSk0qfBRliTknmHu_jHbHxa4U0Wou0gxhPvTeZ_cM/edit#gid=0
02:08:29	Matt Nieto:	¿la matriz u sería como la evolución espacial?
02:16:28	Matt Nieto:	¿porqué la amplitud se hace menos?
02:20:29	Matt Nieto:	no se ve
02:23:12	Matt Nieto:	¿mide la distancia desde el punto a la distancia del origen , cierto?
02:23:41	Matt Nieto:	ya entendí, muchas gracias
02:43:30	Christian Barreto:	BAJEN LOS DATOS DEL DIA 09/08
02:43:45	Christian Barreto:	shape files + excel
02:48:18	Christian Barreto:	#####################################
02:49:41	Christian Barreto:	import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

file="datos_puno_eofs.xlsx"
data=pd.read_excel(file,"DATOS")
location=pd.read_excel(file,"UBICACIONES")
data
02:51:30	Christian Barreto:	data_idx=data.set_index("FECHA")

data_ene=data_idx.loc[data_idx.index.month==1]
data_ene_c=data_ene.loc['1981-01-31':'2010-01-31']
data_ene_c_anomstd=(data_ene_c-data_ene_c.mean())/data_ene_c.std()
data_ene_c_anomstd
02:53:11	Christian Barreto:	data_ene_c_anomstd.plot()
02:53:15	Christian Barreto:	data_ene_c_anomstd["St_114035"].plot()
02:54:31	Christian Barreto:	serite_totest=data_ene_c_anomstd["St_114035"].values
02:55:59	Christian Barreto:	data_ene_c_anomstd=data_ene_c_anomstd.T
02:56:02	Christian Barreto:	data_ene_c_anomstd
02:56:20	Christian Barreto:	matriz_F=data_ene_c_anomstd.to_numpy()
02:57:13	Christian Barreto:	# matris de covarianza, es corrc on anom std, RFF
RFF = (matriz_F @ matriz_F.T) /29 # RFF
print(RFF.shape)
print("------------------")
RFF
02:59:02	Christian Barreto:	# 2) eigenproblem
e_values, e_vectors = np.linalg.eig(RFF)

print("Eigenvalues :\n", e_values)
print("---------------")
print("Eigenvectors:\n", e_vectors)
02:59:26	Christian Barreto:	isort = np.argsort(e_values)[::-1]
e_values = e_values[isort]
e_vectors = e_vectors[:, isort] # re-arrange the columns
print("Eigenvalues :\n", e_values)
print("---------------")
print("Eigenvectors:\n", e_vectors)
03:06:00	Christian Barreto:	A = e_vectors.T @ matriz_F
print(A.shape)
A
03:11:34	Christian Barreto:	print("varianzas:\n", np.var(A,axis=1,ddof=1)) # N-1/ N-ddof
print("eigenvalues:\n", e_values)
print("% de explicacion :\n", 100*e_values/np.sum(e_values))
03:11:59	Christian Barreto:	# plotear componentes principales, evolucion temporal de cada modo
plt.figure()
plt.plot(A.T)
plt.legend(["0", "1", "2", "3","4","5","6","7","8","9","10","11","12","13"])
plt.title("PCs / EOF temporal evolu")
03:15:15	Christian Barreto:	dat_from_eofs = e_vectors @ A
03:17:44	Christian Barreto:	plt.plot(dat_from_eofs[1,:],"r--",label="reconstruido")
plt.plot(serite_totest,label="original",alpha=0.5)
03:19:32	Christian Barreto:	k=3
dat_from_4eofs = e_vectors[:,:k] @ A[:k,:]
plt.plot(np.arange(1981,2011),dat_from_4eofs[1,:],"r--",label="reconstruido")
plt.plot(np.arange(1981,2011),serite_totest,label="original",alpha=0.5)
03:23:17	Christian Barreto:	##
03:23:20	Christian Barreto:	# SVD
03:23:21	Christian Barreto:	u, d, vt = np.linalg.svd(matriz_F, full_matrices=False)

print("Eigenvalues from Covarianza Matrix:\n", e_values)
print("Eigenvalues from SVD:\n", (d ** 2) / 29)  # Note the power of 2. and (N-1)

vecs_svd = u
print("\nEigenvectors from Covarianza Matrix (columns):\n", e_vectors)
print("Eigenvectors from SVD (columns):\n", vecs_svd)
03:24:25	Christian Barreto:	plt.figure()
plt.plot(vt[0,:],"r--",label="PC1 DEL SVD")
plt.plot(A[0,:],label="PC1 DE MATRIX DE COV",alpha=0.5)
plt.legend()
03:24:59	Christian Barreto:	PCsvd = np.dot(np.diag(d), vt)
plt.figure()
plt.plot(PCsvd[0,:]*-1,"r--",label="PC1 DEL SVD")
plt.plot(A[0,:],label="PC1 DE MATRIX DE COV",alpha=0.5)
plt.legend()
03:25:35	Christian Barreto:	for i in range(0,13):
    plt.plot(PCsvd[i,:],label="PC"+str(i+1))
plt.legend()
03:27:21	Christian Barreto:	plt.plot(PCsvd[0,:]/np.sqrt(e_values[0]),"r--",label="PC1 DEL SVD normalizado")
plt.plot(PCsvd[12,:]/np.sqrt(e_values[12]),"b--",label="PC13 DEL SVD normalizado")
plt.legend()
03:28:50	Christian Barreto:	plt.plot(e_values*100/np.sum(e_values),"o-")
03:29:26	Christian Barreto:	# la regla de North et al Thumb
var_errors=e_values*np.sqrt(2/29) # grados de libertad N*
#plt.plot(e_values*100/np.sum(e_values),"o-")
plt.errorbar(np.arange(1,14), e_values*100/np.sum(e_values), yerr=var_errors*100/np.sum(e_values), fmt='xb');
03:30:46	Christian Barreto:	# la regla de North et al Thumb
var_errors=e_values*np.sqrt(2/29) # grados de libertad N*
#plt.plot(e_values*100/np.sum(e_values),"o-")
plt.errorbar(np.arange(1,14), e_values*100/np.sum(e_values), yerr=var_errors*100/np.sum(e_values), fmt='xb');
03:31:48	Christian Barreto:	# correlation PCs vs EOFs

matriz_F_corr=np.zeros([PCsvd.shape[0],PCsvd.shape[0]],dtype=float)

for i in np.arange(0,PCsvd.shape[0]):
    for j in np.arange(0,matriz_F.shape[0]):
        pcinfo=PCsvd[i,:]
        finfo=matriz_F[j,:]
        
        print(np.corrcoef(pcinfo,finfo)[0 , 1])
        
        matriz_F_corr[j,i]=np.corrcoef(pcinfo,finfo)[0 , 1]
    
    print("------------------------")
03:34:03	Christian Barreto:	# mapa de PCs con corr F

import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
import matplotlib.ticker as mticker
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter


fname = 'INEI_LIMITE_DEPARTAMENTAL_GEOGPSPERU_JUANSUYO_931381206.shp'
shape_feature = ShapelyFeature(Reader(fname).geometries(),
                              ccrs.PlateCarree(), facecolor='none',edgecolor='k',linewidth=0.5 )

plain_crs = ccrs.PlateCarree()
puno_extent = [-71.25, -68.7, -17.4,-12.9]
to_test = [(plain_crs, puno_extent)]

projection=plain_crs
extent=puno_extent
03:34:07	Christian Barreto:	fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(1, 1, 1, projection=projection)
ax.add_feature(shape_feature)

ax.set_extent(extent, crs=plain_crs)
03:34:12	Christian Barreto:	ax.set_xticks([-71,-70,-69], crs=ccrs.PlateCarree())
ax.set_yticks([-17.5,-17,-16.5,-16,-15.5,-15,-14.5,-14,-13.5,-13], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,
                                   number_format='.1f')
lat_formatter = LatitudeFormatter(number_format='.1f')
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
03:34:19	Christian Barreto:	x=location.LON.values.tolist()
y=location.LAT.values.tolist()
03:34:23	Christian Barreto:	plt.scatter(x, y, c=list(matriz_F_corr[:,0]),s=200,cmap='bwr',vmin=-1, vmax=1)
03:34:28	Christian Barreto:	plt.title("EOF-1 with corr")
plt.colorbar()
03:37:57	Christian Barreto:	####################
03:37:58	Christian Barreto:	# mapa de PCs con corr F

import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
import matplotlib.ticker as mticker
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter


fname = 'INEI_LIMITE_DEPARTAMENTAL_GEOGPSPERU_JUANSUYO_931381206.shp'
shape_feature = ShapelyFeature(Reader(fname).geometries(),
                              ccrs.PlateCarree(), facecolor='none',edgecolor='k',linewidth=0.5 )

plain_crs = ccrs.PlateCarree()
puno_extent = [-71.25, -68.7, -17.4,-12.9]
to_test = [(plain_crs, puno_extent)]
03:38:03	Christian Barreto:	projection=plain_crs
extent=puno_extent

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(1, 1, 1, projection=projection)
ax.add_feature(shape_feature)

ax.set_extent(extent, crs=plain_crs)
03:38:07	Christian Barreto:	ax.set_xticks([-71,-70,-69], crs=ccrs.PlateCarree())
ax.set_yticks([-17.5,-17,-16.5,-16,-15.5,-15,-14.5,-14,-13.5,-13], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,
                                   number_format='.1f')
lat_formatter = LatitudeFormatter(number_format='.1f')
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
03:38:12	Christian Barreto:	#---------------------------
# leemos la data de longitud y latitud
#---------------------------
#data_st=results_alltrends[results_alltrends['Tendencia'] == 'Tendencia negativa no significativa' ]
x=location.LON.values.tolist()
y=location.LAT.values.tolist()

#plt.scatter(x, y)
plt.scatter(x, y, c=list(matriz_F_corr[:,1]),s=200,cmap='bwr',vmin=-1, vmax=1)
            #,zorder=2,marker="v",color="none", edgecolor="blue",label='Tend. negativa no sig.')
    
plt.title("EOF-2 with corr")
plt.colorbar()
03:39:56	Christian Barreto:	pip install eofs
03:43:17	Christian Barreto:	##########################
03:43:21	Christian Barreto:	#ahora con sst
03:43:42	Christian Barreto:	import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

dset = xr.open_dataset("sst.mnmean.nc")
dset
03:43:55	Christian Barreto:	sst=dset.sst
sst_trop=sst.sel(lat=slice(30, -30))
03:44:45	Christian Barreto:	month_idxs=sst_trop.groupby('time.month').groups
jan_idxs=month_idxs[1]
sst_trop_jan=sst_trop.isel(time=jan_idxs)
sst_trop_jan=sst_trop_jan.sel(time=slice('1981-01-01', '2010-01-01'))
# seleccionamos solo el rango climatico
sst_trop_jan_c=sst_trop_jan.sel(time=slice('1981-01-01', '2010-01-01')).mean(dim='time',skipna=True)
sst_trop_jan_std=sst_trop_jan.sel(time=slice('1981-01-01', '2010-01-01')).std(dim='time',skipna=True)
03:45:36	Christian Barreto:	sst_trop_jan_anomstd=(sst_trop_jan-sst_trop_jan_c)/sst_trop_jan_std #anomalias std
03:49:31	Christian Barreto:	from eofs.xarray import Eof
03:49:34	Christian Barreto:	solver = Eof(sst_trop_jan_anomstd,center=False)
03:49:37	Christian Barreto:	sst_eofs = solver.eofsAsCorrelation()
sst_pcs = solver.pcs(npcs=3, pcscaling=1)
03:49:41	Christian Barreto:	sst_eofs[0].plot(robust=True)
plt.title("modo 1")
03:51:19	Christian Barreto:	sst_pcs[:,0].plot(marker="o") # modo 1
plt.title("modo 1")

plt.grid()
03:52:59	Christian Barreto:	###############################
sst_eofs[1].plot(robust=True)
plt.title("modo2")
03:53:03	Christian Barreto:	sst_pcs[:,1].plot(marker="o") # modo 1
plt.title("modo2")
plt.grid()
03:54:19	Christian Barreto:	## modo 3
sst_eofs[2].plot(robust=True)
plt.title("modo3")
03:54:23	Christian Barreto:	sst_pcs[:,2].plot(marker="o") # modo 1
plt.title("modo3")
plt.grid()
03:54:32	Christian Barreto:	#https://ajdawson.github.io/eofs/latest/api/eofs.standard.html
03:55:06	Christian Barreto:	## varianzas por cada PCs
sum_ev=solver.totalAnomalyVariance()
var_pcs=solver.varianceFraction()
plt.plot(var_pcs*100,marker="*")
03:55:41	Christian Barreto:	# north et al rule thumb
pcs_var_errors=solver.northTest()
errorbar=(pcs_var_errors*100/sum_ev).values
plt.errorbar(np.arange(1,31), var_pcs*100, yerr=errorbar, fmt='xb');
03:56:18	Christian Barreto:	# north et al rule thumb
pcs_var_errors=solver.northTest()
errorbar=(pcs_var_errors*100/sum_ev).values
plt.errorbar(np.arange(1,31), var_pcs*100, yerr=errorbar, fmt='xb');
04:04:29	Christian Davila:	No tengo inconveniente con que sea viernes
04:04:32	Angela Peña:	JUEVES
04:04:32	Stefany Amado Menauth DZ6:	Viernes
04:04:35	Rosamaría Pérez Bellido:	viernes estaría bien
04:04:39	Angela Peña:	O VIERNES
04:04:41	Jorge Poma  Dz11:	igual jueves o viernes :)
04:04:44	IMELDA VALENTINA ALIAGA GUERREROS:	VIERNES
04:04:47	Gabriel Leví Caro-Sánchez Gago:	ok viernes
04:04:49	LEONARDO MESIAS CHAVEZ LOPEZ:	viernes
04:04:52	Lia Nicolle:	claro viernes
04:04:54	sandra rivadeneria:	viernes porf
04:05:11	Angela Peña:	SI