祁连山综合观测网：黑河流域地表过程综合观测网（混合林站自动气象站-2018）

import pandas as pd
df=pd.read_excel(r'C:\Users\mengx\Desktop\土壤水分的频谱分析\祁连山综合观测网：黑河流域地表过程综合观测网（混合林站自动气象站-2018）\2018年黑河流域地表过程综合观测网混合林站AWS.xlsx')

将所有包含-6999异常值的数据行都删除

df=df[~df.isin([-6999])].dropna()
df.describe()

	WS_28m	WD_28m	Ta_28m	RH_28m	Rain	Press	IRT_1	IRT_2	PAR_down	PAR_up	...	Ts_2cm	Ts_4cm	Ts_10cm	Ts_20cm	Ts_40cm	Ts_60cm	Ts_100cm	Ts_160cm	Ts_200cm	Ts_240cm
count	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	...	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000	48710.000000
mean	4.077554	175.662670	11.461099	29.816694	0.001866	916.709485	13.927786	14.165521	426.747442	32.581842	...	10.470963	11.022871	8.774166	10.764586	10.534167	10.350921	10.117353	10.088647	10.135128	13.152512
std	2.025049	95.164170	14.144133	18.637659	0.088265	9.250470	15.527007	15.768607	585.469889	45.463257	...	13.383667	12.179141	9.838000	10.818108	10.246014	9.131619	7.594671	5.522301	4.868285	4.268635
min	0.000000	0.000000	-26.491390	0.550920	0.000000	894.720520	-25.433640	-24.563670	0.000000	0.000000	...	-19.066890	-13.060090	-19.140970	-7.747810	-17.799430	-1.805150	0.301300	2.470000	3.355000	7.149000
25%	2.913832	76.387883	1.011122	15.444168	0.000000	909.384980	2.474515	2.608615	0.000000	0.000000	...	0.139390	1.109818	0.102010	0.440200	0.413803	-0.188868	2.079925	4.588000	5.312000	9.640000
50%	3.803475	204.708120	13.568130	25.073510	0.000000	916.599215	15.223605	15.288935	26.896695	1.065090	...	11.535135	11.834830	9.176530	11.085925	11.208225	10.599165	10.092180	9.870000	9.850000	11.890000
75%	4.955337	258.739375	22.890038	39.274338	0.000000	923.325688	25.739428	26.046557	800.482710	61.899407	...	22.175790	22.660810	18.135233	21.464008	20.839128	19.362520	17.381982	15.320000	15.200000	17.240000
max	16.303840	359.204320	37.903010	100.000000	13.700000	948.232790	51.118450	59.240230	2451.167240	275.145080	...	37.444560	33.805460	29.103050	33.129830	29.759500	27.538020	26.394470	19.910000	17.180000	21.910000

8 rows × 39 columns

#取出全部的含水率数据
ms=df.iloc[:,19:29]

import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] # 解决汉字显示为□指定默认字体为黑体。
plt.rcParams['axes.unicode_minus'] = False
#不同深度的含水率变化
ms.plot()
plt.legend(loc=2, bbox_to_anchor=(1.05,0.7),borderaxespad = 0.)
plt.title(r'各深度的土壤水分')
plt.grid()

import matplotlib.pyplot as plt
# plt.plot(df["Ms_160cm"],'-k')
# plt.plot(df["Ms_4cm"],'-r')
# plt.title('4cm和160cm深度的土壤水分')
plt.plot(df["Rain"],'--g')
plt.grid()
plt.show()

#看一看谱密度
import pandas as pd
df=pd.read_excel(r'C:\Users\mengx\Desktop\2018年黑河流域地表过程综合观测网混合林站AWS.xlsx')
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
# 删除-6999
df=df[~df.isin([-6999])].dropna()
p,f=plt.psd(df["Ms_4cm"], NFFT=10000, Fs=2, window=mlab.window_none, pad_to=None,return_line=False)
plt.title('4cm的功率谱密度')

Text(0.5, 1.0, '4cm的功率谱密度')

#看一看谱密度
import pandas as pd
df=pd.read_excel(r'C:\Users\mengx\Desktop\土壤水分的频谱分析\2018年黑河流域地表过程综合观测网混合林站AWS.xlsx')
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['STSong']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False  # 用来正常显示负号
plt.rcParams['font.size']=14#设置全局字体
import matplotlib.mlab as mlab
# 删除-6999
df=df[~df.isin([-6999])].dropna()
p,f=plt.psd(df["Rain"], NFFT=5096, Fs=1, window=mlab.window_none, pad_to=None,return_line=False)
plt.title('降雨的功率谱密度')
fig2=plt.figure(2)
plt.loglog(f,p)
plt.grid()

df.shape

(48710, 40)

df=pd.concat([pd.DataFrame(f),pd.DataFrame(p)],axis=1)
df=df[~df.isin([0])].dropna()

import numpy as np
plt.loglog(df.iloc[:,0],df.iloc[:,1])
plt.grid()
from scipy import optimize
def f_1(x, A, B):
    return A * x + B
A1, B1 = optimize.curve_fit(f_1, np.log10(df.iloc[1000:,0]),np.log10(df.iloc[1000:,1]))[0]
def f_2(x, A, B):
    return A * x + B
A2, B2 = optimize.curve_fit(f_2, np.log10(df.iloc[:400,0]),np.log10(df.iloc[:400,1]))[0]
print(A1,B1)
print(A2,B2)
x1=np.arange(np.log10(df.iloc[:,0]).min(),np.log10(df.iloc[:,0]).max(),0.02)
# print(x1)
y1 =f_1(x1,A1,B1)
y2 =f_2(x1,A2,B2)
#绘图
plt.loglog(10**x1,10**(y1))
plt.loglog(10**x1,10**(y2))
plt.text(10e-4,10e-3,r'$\beta 1=$'+str(A1)[1:7])
plt.text(10e-4,10e-1,r'$\beta 2=$'+str(A2)[1:7])
plt.show()

-0.9056930734066293 -2.49921212944677
-0.13419953778885224 -1.6300723085417057

黑河上游生态水文传感器网络逐时土壤水分观测数据集（2013-2017）

import pandas as pd
four_year=pd.read_csv(r'C:\Users\mengx\Desktop\土壤水分的频谱分析\2013-2018_all.csv').dropna(axis=0)

import matplotlib.pyplot as plt
plt.figure(figsize=(15,3))
plt.plot(four_year.iloc[:,2])
plt.show()

import matplotlib.mlab as mlab
p,f=plt.psd(four_year.iloc[:,2], NFFT=20000, Fs=1, window=mlab.window_none, pad_to=None,return_line=False)

删除0元素所在的行

df=pd.concat([pd.DataFrame(f),pd.DataFrame(p)],axis=1)
df=df[~df.isin([0])].dropna()

拟合对数坐标下的斜率

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['STSong']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False  # 用来正常显示负号
plt.rcParams['font.size']=14#设置全局字体
plt.loglog(df.iloc[:,0],df.iloc[:,1])
plt.grid()
from scipy import optimize
def f_1(x, A, B):
    return A * x + B
A1, B1 = optimize.curve_fit(f_1, np.log10(df.iloc[:,0]),np.log10(df.iloc[:,1]))[0]
print(A1,B1)
x1=np.arange(np.log10(df.iloc[:,0]).min(),np.log10(df.iloc[:,0]).max(),0.02)
# print(x1)
y1 =f_1(x1,A1,B1)
#绘图
plt.loglog(10**x1,10**(y1))
plt.text(10e-4,10e-4,r'$\beta=$'+str(A1)[:7])
plt.show()

-2.065160490229184 -5.742930857784623