import pandas as pdfrom standard_precip import spi
df_rainfall = pd.read_csv('/home/mw/project/wichita_rain.csv')df_rainfall.head()
spi_rain = spi.SPI()spi_monthly = spi_rain.calculate(df_rainfall, 'date', 'precip', freq="M", scale=1,
fit_type="lmom", dist_type="gam")spi_monthly.head()
spi_3monthly = spi_rain.calculate(df_rainfall, 'date', 'precip', freq="M", scale=3,
fit_type="mle", dist_type="gam")
spi_3monthly.head()
import matplotlib.pyplot as pltfig = plt.figure(figsize = (10,8))# 绘图# 获取日期的索引并选取特定间隔的日期作为标签date_indices = spi_3monthly.index[::120] # freq 是你的频率，例如每隔7天选取一个日期selected_dates = spi_3monthly.loc[date_indices, 'date']plt.plot(spi_3monthly.date, # x轴数据spi_3monthly.precip_scale_3_calculated_index, # y轴数据linestyle = '-', # 折线类型linewidth = 2, # 折线宽度color = 'steelblue', # 折线颜色marker = 'o', # 点的形状markersize = 6, # 点的大小markeredgecolor='black', # 点的边框色markerfacecolor='brown') # 点的填充色# 添加标题和坐标轴标签plt.title('SPI3')plt.xlabel('日期')plt.ylabel('SPI')# 剔除图框上边界和右边界的刻度plt.tick_params(top = 'off', right = 'off')# 显示图形plt.show()
from gma import climetData = pd.read_csv('/home/mw/project/wichita_rain.csv')PRE = Data['precip'].values# 分别计算1个月、3个月、6个月、12个月、24个月、60个月尺度的 SPI 指数SPI1 = climet.Index.SPI(PRE)SPI3 = climet.Index.SPI(PRE, Scale = 3)SPI6 = climet.Index.SPI(PRE, Scale = 6)SPI12 = climet.Index.SPI(PRE, Scale = 12)SPI24 = climet.Index.SPI(PRE, Scale = 24)SPI60 = climet.Index.SPI(PRE, Scale = 60)
data = {'date': Data['date'],'SPI1': SPI1,'SPI3': SPI3,'SPI6': SPI6,'SPI12': SPI12,'SPI24': SPI24,'SPI60': SPI60,}# 将数据转换为DataFramedf = pd.DataFrame(data)df
fig = plt.figure(figsize = (12,8))# 绘图# 获取日期的索引并选取特定间隔的日期作为标签date_indices = df.index[::120] # freq 是你的频率，例如每隔7天选取一个日期selected_dates = df.loc[date_indices, 'date']plt.plot(df.date, # x轴数据df.SPI3, # y轴数据linestyle = '-', # 折线类型linewidth = 2, # 折线宽度color = 'steelblue', # 折线颜色marker = 'o', # 点的形状markersize = 6, # 点的大小markeredgecolor='black', # 点的边框色markerfacecolor='brown') # 点的填充色# 添加标题和坐标轴标签plt.title('SPI3')plt.xlabel('日期')plt.ylabel('SPI')plt.xticks(date_indices, selected_dates, rotation=45)# 剔除图框上边界和右边界的刻度plt.tick_params(top = 'off', right = 'off')# 显示图形plt.show()
