00:44:35	Christian Barreto:	%reset -sf

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

file="datos_reg_lineal_esta.xlsx"
data=pd.read_excel(file)

data
00:45:01	Christian Barreto:	y=data["chuqui_jja"]
x=data["chuqui_mam"]
plt.scatter(x,y)
00:45:49	Christian Barreto:	plt.plot(x)
00:46:02	Christian Barreto:	plt.plot(y)
00:47:59	Christian Barreto:	#create pandas para acumular datos

df= pd.DataFrame({"x":x,"y": y})

#calculamos la media de x e y

xmean= np.mean(x)
ymean= np.mean(y)

#calculamos los terminos que necesitamos en el numerador y denominador de beta
df["xycov"] = (df["x"] - xmean) * (df["y"] - ymean)
df["xvar"] = (df["x"] - xmean)**2

#calculamos a y b
b=df["xycov"].sum()/ df["xvar"].sum()
a= ymean-(b*xmean)
print("a:",a)
print("b:",b)
00:48:53	Christian Barreto:	ypred = a  + b *x
00:52:00	Christian Barreto:	# plotear la linea
plt.figure()
plt.plot(x, ypred) # regresion line
plt.plot(x,y,"ro") # scatter plot
plt.title("r lineal")

plt.xlabel("chuqui_mam")
plt.ylabel("chuqui_jul")
00:55:02	Christian Barreto:	# usando sklearn 

from sklearn.linear_model import LinearRegression

lm = LinearRegression()
model = lm.fit(x.to_numpy().reshape((-1,1)),y.to_numpy())
00:56:35	Christian Barreto:	print(f"alpha = {model.intercept_}")
print(f"beta = {model.coef_}")
01:01:55	Christian Barreto:	plt.plot(np.arange(1971,2019),y,label="yobs")
plt.plot(np.arange(1971,2019),ypred,label="ypred")
plt.legend()
01:14:05	Christian Barreto:	print("std pred:", np.std(ypred))
print("std obs:", np.std(y))
corr_matrix=np.corrcoef(ypred,y)
correlation_xy =corr_matrix[0,1]
r_squared = correlation_xy**2
print("r**2:", r_squared)
01:31:39	Christian Barreto:	##############################
01:31:42	Christian Barreto:	tednencias
01:31:46	Christian Barreto:	%reset -sf
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
file="DATA_TENDENCIAS.xlsx"
data=pd.read_excel(file)
data
01:35:04	Christian Barreto:	data_idx=data.set_index("DATE")
data_ene=data_idx.loc[data_idx.index.month==1]
data_ene
01:43:25	Christian Barreto:	data_ene_np=data_ene.TARMA_TX.values
01:44:40	Christian Barreto:	from sklearn.linear_model import LinearRegression

lm = LinearRegression()

x=np.arange(1964,2019)
y=data_ene_np

model = lm.fit(x.reshape((-1,1)),y)

b=model.coef_[0]
a=model.intercept_
linea= a + b*x

plt.plot(x,y)
plt.plot(x,linea,"r--")
01:58:45	Christian Barreto:	# evolucion temporal de tendencias

# del 70 al 00
x1=np.arange(1970,2001)
y1=data_ene.TARMA_TX.loc["1970-01-31":"2000-01-31"].values

model =lm.fit(x1.reshape((-1,1)),y1)
b1=model.coef_[0]
a1=model.intercept_
linea1 = a1 + b1 *x1

# del 2000 al 2018

x3=np.arange(2000,2019)
y3=data_ene.TARMA_TX.loc["2000-01-31":"2018-01-31"].values
model =lm.fit(x3.reshape((-1,1)),y3)
b3=model.coef_[0]
a3=model.intercept_
linea3 = a3 + b3 *x3
01:58:48	Christian Barreto:	#plots
plt.plot(x,y)
plt.plot(x,linea,"r--",label="pend.todo= "+str(np.round(b,3)))
plt.plot(x1,linea1,"k--",label="pend.70-00= "+str(np.round(b1,3)))
plt.plot(x3,linea3,"k--",label="pend.00-18= "+str(np.round(b3,3)))
plt.legend()
02:55:35	Christian Barreto:	####################################
02:55:37	Christian Barreto:	###################################
02:55:38	Christian Barreto:	import pymannkendall as mk

result=mk.original_test(y,alpha=0.05)

print("slope=",result.slope)
print("pvalor=",result.p)
02:56:41	Christian Davila:	pip install pymannkendall
02:56:56	Christian Barreto:	conda install -c conda-forge pymannkendall
03:04:38	Christian Davila:	f"pend= {np.round(slp,2)} {np.round(pv,6)}"

Puedes ponerlo de esta forma para evitar conflicto al concatenar
03:06:39	Christian Barreto:	r = mk.original_test(y,alpha=0.05)
pv= r.p
itr= r.intercept
slp = r.slope

x=np.arange(1,56) #1 a 55
lineakd = itr + slp*x

plt.plot(x,y,marker="o")
plt.plot(x,lineakd,label="pend="+str(np.round(slp,3))+" y p-val:"+str(np.round(pv,6)))
plt.legend()

pos=np.arange(1,56,10)
plt.grid()
labels=np.arange(1964,2019,10)
plt.xticks(pos,labels)
03:13:17	Christian Barreto:	# 90 - 2018
y2=data_ene.TARMA_TX.loc['1991-01-31':'2018-01-31'].values
r2 = mk.original_test(y2,alpha=0.05)
pv2 = r2.p
print(pv2)
itr2 = r2.intercept
print(itr2)
slp2= r.slope

x2=np.arange(28,55) #1a 55
lineakd2 = itr2 + slp2 * x2


#plots
plt.plot(x,y,marker="o")
plt.plot(x2,lineakd2,label="pend="+str(np.round(slp2,3))+" y p-val:"+str(round(pv2,6)))
plt.legend()
03:19:50	ALANIA SUMARAN JOEL YOEL DZ02:	precip.mon.mean
03:19:52	ALANIA SUMARAN JOEL YOEL DZ02:	este??
03:27:33	Christian Barreto:	########################
03:27:37	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
03:27:42	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
jan_idxs=month_idxs[1]

#extramos solo los eneros

data=data_c.isel(time=jan_idxs)
data
03:33:15	Christian Barreto:	import pymannkendall as mk

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

for ni in range(0,data.shape[2]):
    for nj in range(0,data.shape[1]):
        
        info=data.isel(lat=nj,lon=ni).values
        
        array_sum=np.sum(info)
        array_has_nan = np.isnan(array_sum)
        
        if array_has_nan == True:
            
            slopes[nj,ni] =np.nan
            pvalues[nj,ni] =np.nan
        
        else:
            
            result = mk.original_test(info,alpha=0.05)
            
            slopes[nj,ni]=result.slope
            pvalues[nj,ni] = result.p
03:35:59	Christian Barreto:	data_set = xr.Dataset( coords={'lon': ([ 'lon'], data.lon),
                                             'lat': (['lat',], data.lat)})

data_set["slope"] = (["lat", "lon"], slopes)
data_set["pvalues"] = (["lat", "lon"], pvalues)

data_set.to_netcdf("trend_81-10.nc", mode="w")
03:37:25	Christian Barreto:	#ploteamos tendencias en pp en enero
dset_2=xr.open_dataset("trend_81-10.nc")

tendencia=dset_2.slope
pv=dset_2.pvalues

tendencia.plot(robust=True)
03:39:48	Christian Barreto:	##################
03:39:50	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)
03:39:55	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)

#adding grillas
03:40:02	Christian Barreto:	ax.gridlines(draw_labels=False, xlocs=xticks, ylocs=yticks)

#plot data
levels=[-2.5,-2,-1,0,1,2,3]

####
tendencia.plot(robust=True)

plt.title("tendencia de pp 81-10")

ax.set_aspect("auto",adjustable=None)
03:59:12	Christian Barreto:	###############################################
03:59: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
03:59:24	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
jan_idxs=data_c.isel(time=jan_idxs)
jan_idxs
04:02:06	Christian Barreto:	############
04:02:07	Christian Barreto:	https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php
04:02:12	Christian Barreto:	#---------------
04:02:20	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]
04:03:53	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
04:04:53	Christian Barreto:	nino_idx=[]
nina_idx=[]
normal_idx=[]

year_idxs=data_jan.groupby("time.year").groups
04:06:04	Christian Barreto:	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])
04:07:55	Christian Barreto:	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)
04:11:03	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)
04:11:07	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)

#adding grillas
04:11:11	Christian Barreto:	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)
04:14:31	Stefany Amado Menauth DZ6:	Gracias