作业要求：

目的：学习用地面站点实测降雨数据校正卫星遥感降雨数据。
原始数据：
2007-2013年的TRMM卫星3B43产品（2013年12个月的月降雨数据，覆盖范围精度-180_{180、纬度-50}50，0.25x0.25度分辨率，像素值单位：mm/hour）；
华北8个省市的气象站2007-2013年的年降雨量文件，以及气象站的SHAPE文件；
数据处理要求：
（1）实验区1（华北8个省市）
对覆盖范围（经度：111.25-121.25,纬度：36-43）内的卫星年降雨数据，用地面数据进行校正，并存为TIF文件。
校正的卫星年降雨数据可以是一年的，如：2010年或2013年，也可以是2007-2013年时间序列数据。

详细过程

选取2009-2013年的年降水数据，对京津冀及周边地区的大概范围（经度：111.25-121.25，纬度：36-43）进行影像降水数据的校正。前提假设站点获得的降水数据为真值。另外，增加了全国dem数据，提高校正精度。总体思路如下：
（1） 读取2009-2013年，研究区内各站点的年降水量、经度、纬度、DEM、站点所处位置对应的遥感数据TRMM3B43观测得到的降水量。
（2） 利用2009-2012年的数据，分析站点的年降水量和经度、纬度、DEM、遥感数据TRMM3B43观测得到的降水量之间的相关性，同时，建立多元回归函数。
（3） 利用2013年的遥感数据TRMM3B43观测得到的降水量，经度、纬度、DEM计算校正后的降水量，该降水量与2013年站点降水数据比较，进行误差分析，评估校正效果。

图1 总体技术路线

import os
from osgeo import gdal
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LinearRegressionif __name__ == '__main__':# 第一步：prec_p = r"D:\prec_corr_2021\input\prec"surf_p = r"D:\prec_corr_2021\input\SURF_CLI_CHN_MUL_YER_STATION.xls"S_prec200912_ch, S_prec2013_ch,S_cor = sta_pre(prec_p, surf_p)#第二步：TRMM_p = "D:\\prec_corr_2021\\input\\RS\\"T_prec200912_ch, T_prec2013_ch= TRMM_pre(TRMM_p,S_cor)#需要遥感提取数据和站点一一对应ST_prec200912,ST_prec2013=[],[]for T_name, T_time, T_prec in T_prec200912_ch:for s_time, s_name, s_prec, s_x, s_y, s_z in S_prec200912_ch:if T_name == s_name and T_time==s_time:ST_prec200912.append((s_time, s_name, s_prec, T_prec,s_x, s_y, s_z))for s_time, s_name, s_prec, s_x, s_y, s_z in S_prec2013_ch:for T_name, T_time, T_prec in T_prec2013_ch:if T_name == s_name and T_time == s_time:ST_prec2013.append((s_time, s_name, s_prec, T_prec, s_x, s_y, s_z))#第三步：a,b = f_cor(ST_prec200912)print('多元回归函数为：')print('y=',a,'+',b[0],'*','T_prec','+',b[1],'*','lon',b[2],'*','lat','+',b[3],'*','dem')#第四步：MAE1, R1,RMSE1,MAE2, R2,RMSE2 = S_val(ST_prec2013,a,b)print('2013年站点数据与校正后遥感数据的误差分析：')print('校正前：平均绝对误差',MAE1,'均方根误差',RMSE1,'相关系数',R1)print('校正后：平均绝对误差',MAE2,'均方根误差',RMSE2,'相关系数',R2)#第五步：TRMM_p = "D:\\prec_corr_2021\\input\\RS\\"DEM_p = "D:\\prec_corr_2021\\input"T_val(DEM_p, TRMM_p,a,b)print("最终结果查看：D:\prec_corr_2021\code")print('-----Complete！-----')

第一个函数主要来实现TXT数据的批量读取以及excel数据的读取，可以有更简单的读取方式，可以自己改进。

#------------------1、read the station data by --------------------
'''
说明：读取气象站降水数据及预处理
输入：降雨数据所在地址（D:\prec_corr_2021\input\prec），气象站综合信息表（SURF_CLI_CHN_MUL_YER_STATION）
输出：站点数组（站点S_name，经度S_x，纬度S_y，海拔S_z）；2009-2012站点降水数据（时间S_time，站点S_name，降水S_prec）;2013站点降水（同上）
'''
def sta_pre(prec_path, surf_path):S_cor, S_prec200912, S_prec2013 = [], [], []if not os.path.exists(prec_path) or not os.path.exists(prec_path):print("prec_path or surf_path is not exist.")return [], [], []if not os.path.isdir(prec_path):print("prec_path is not folder.")return [], [], []# 读取全部的站点数据try:df = pd.read_excel(surf_path,sheet_name="Sheet1",usecols=[0, 11, 12, 13])  # 指定读取第1/12/13/14列except Exception:print("read excel failed")return [], [], []# print(df)# 获取{站点S_name:(站点S_name，经度S_x，纬度S_y，海拔S_z)}df = np.array(df)dic_df = {}for item in df:s_name_ex, s_x, s_y, s_z = item.tolist()# 只有站点名称是浮点数字才认为是正确的if isinstance(s_name_ex, float):s_name_ex = str(int(s_name_ex))dic_df[s_name_ex] = (s_name_ex, s_x, s_y, s_z)# print(len(dic_df))# 读取站点降水数据file_names = os.listdir(prec_path)file_content = []for filename in file_names:file_path = os.path.join(prec_path, filename)with open(file_path, "r") as f:file_content = file_content + f.readlines()def get_info(l_line, years):stime, sname, sprec = l_linestime = stime.strip()if stime in years:try:# 做一个判断，剔除降水值为负的值i_prec = int(sprec.strip())if i_prec > 0:return stime.strip(), sname.strip(), i_precexcept Exception:return "", "", 0return "", "", 0# print(file_content)for line in file_content:list_line = line.strip().split(" ")if len(list_line) == 3:s_time, s_name, i_prec = get_info(list_line, ["2009", "2010", "2011", "2012"])if s_time != "" or s_name != "":S_prec200912.append((s_time, s_name, i_prec))s_time, s_name, i_prec = get_info(list_line, ["2013"])if s_time != "" or s_name != "":S_prec2013.append((s_time, s_name, i_prec))# 剔除站点数据中，不在S_prec200912 和 S_prec2013中的数据for s_time, s_name, s_prec in S_prec200912:if s_name in dic_df:S_cor.append(dic_df.get(s_name))for s_time, s_name, s_prec in S_prec2013:if s_name in dic_df:S_cor.append(dic_df.get(s_name))# 去重S_cor = list(set(S_cor))#将经纬度坐标的赋值给站点数据S_prec200912_ch, S_prec2013_ch = [], []for s_name_ex, s_x, s_y, s_z in S_cor:for s_time_1, s_name_1, s_prec_1 in S_prec200912:if s_name_1 == s_name_ex:s_prec_1=s_prec_1*0.1s_prec_1=round(s_prec_1, 2)s_z=round(s_z, 2)s_y=round(s_y, 2)s_z=round(s_z, 2)S_prec200912_ch.append((s_time_1,s_name_1,s_prec_1,s_x, s_y, s_z))for s_time_2, s_name_2, s_prec_2 in S_prec2013:if s_name_2 == s_name_ex:s_prec_2 = s_prec_2 * 0.1s_prec_2=round(s_prec_2, 2)s_z=round(s_z, 2)s_y=round(s_y, 2)s_z=round(s_z, 2)S_prec2013_ch.append((s_time_2, s_name_2, s_prec_2, s_x, s_y, s_z))return S_prec200912_ch, S_prec2013_ch,S_cor

第二个函数主要用来读取遥感数据以及提取站点所在地遥感数据对应的降水量，需要具有一定的GIS基础。但缺少批量读取遥感数据的步骤，可以参考函数1进行改进。

#------------------2、read the TRMM 3B43 data by--------------------
'''
说明：读取遥感卫星降水数据及预处理，获得站点所在遥感数据的值；同时读取对应DEM数据
输入：数据所在地址（D:\prec_corr_2021\input\RS），站点降水数据（时间S_time，站点S_name，降水S_prec）
输出：2009-2012遥感降水数据（T_prec，需和S_prec200912数组一一对应）;2013遥感降水数据（T_prec2013）
'''
def TRMM_pre(TRMM_path,S_cor):if not os.path.exists(TRMM_path) or not os.path.exists(TRMM_path):print("TRMM_path or surf_path is not exist.")return [], [], []if not os.path.isdir(TRMM_path):print("TRMM_path is not folder.")return [], [], []#----------------------①读取TRMM卫星影像数据及DEM数据 - ---------dataset1 = gdal.Open(TRMM_path + '3B43_2009.TIF')dataset2 = gdal.Open(TRMM_path + '3B43_2010.TIF')dataset3 = gdal.Open(TRMM_path + '3B43_2011.TIF')dataset4 = gdal.Open(TRMM_path + '3B43_2012.TIF')dataset5 = gdal.Open(TRMM_path + '3B43_2013.TIF')samples = dataset1.RasterXSizelines = dataset1.RasterYSize# bands = dataset1.RasterCountimg_geotrans = dataset1.GetGeoTransform()# img_proj = dataset1.GetProjection()im_data1 = dataset1.ReadAsArray(0, 0, samples, lines)im_data2 = dataset2.ReadAsArray(0, 0, samples, lines)im_data3 = dataset3.ReadAsArray(0, 0, samples, lines)im_data4 = dataset4.ReadAsArray(0, 0, samples, lines)im_data5 = dataset5.ReadAsArray(0, 0, samples, lines)del dataset1del dataset2del dataset3del dataset4del dataset5# -----------------------②转换站点经纬度坐标到对应TRMM遥感图像上的行列号----------im_loc = [[], []]lon = [i[1] for i in S_cor]lat=[i[2] for i in S_cor]pnt_coordinates = [lon]+[lat]# 第一个列表放行号，存纬度换出来的值，后者为列号，存经度换出来的值for i in range(len(pnt_coordinates[0])):sample = sp.Symbol('sample')line = sp.Symbol('line')a = -pnt_coordinates[0][i] + img_geotrans[0] + sample*img_geotrans[1] + line*img_geotrans[2]b = -pnt_coordinates[1][i] + img_geotrans[3] + sample*img_geotrans[4] + line*img_geotrans[5]answer = sp.solve([a,b],[sample,line])#解二元一次方程组line = int(np.floor(answer[line]))sample = int(np.floor(answer[sample]))im_loc[0] += [line]im_loc[1] += [sample]#----------------------③读取2009_2013年站点对应的TRMM卫星影像数据-------------------------data_TRMM_stationpnt = [[],[],[],[],[]]for i in range(len(im_loc[0])):data_TRMM_stationpnt[0] += [im_data1[im_loc[0][i], im_loc[1][i]]]data_TRMM_stationpnt[1] += [im_data2[im_loc[0][i], im_loc[1][i]]]data_TRMM_stationpnt[2] += [im_data3[im_loc[0][i], im_loc[1][i]]]data_TRMM_stationpnt[3] += [im_data4[im_loc[0][i], im_loc[1][i]]]data_TRMM_stationpnt[4] += [im_data5[im_loc[0][i], im_loc[1][i]]]T_prec200912 = data_TRMM_stationpnt[0:4]T_prec2013 = data_TRMM_stationpnt[4]T_prec200912_ch,T_prec2013_ch=[], []time_y=['2009','2010','2011','2012','2013']for i,(s_name, s_x, s_y, s_z) in enumerate(S_cor):for j,t in enumerate (time_y):if t=='2013':T_prec2013_ch.append((s_name, t, T_prec2013[i]))else:T_prec200912_ch.append((s_name, t, T_prec200912[j][i]))return T_prec200912_ch, T_prec2013_ch

第三个函数主要构建校正函数，这里选择的多元一次函数，可以换成其他方法，如随机森林等。

#------------------3、built model between TRMM and station --------------------
'''
说明：建立2009-2012年，站点数据和TRMM数据T_prec,站点所在经度、纬度、高程的回归关系
输入：2009-2012年站点及遥感降水数据
输出：多元一次回归函数截距及回归系数
'''
def f_cor(S_prec200912):S_prec200912 = np.array(S_prec200912,dtype = float)y =  S_prec200912[:,2]                #实际降水即站点降水x =  S_prec200912.T[np.arange(3,7)]   #影响降水的因子x=x.Tall_data =  S_prec200912.T[np.arange(2,7)]all_data=all_data.Tall_data=pd.DataFrame(all_data,columns=['S_prec','T_prec','lon', 'lat', 'dem'])all_data.boxplot()plt.savefig("01影响因子分析.jpg",dpi=100)# plt.show()plt.clf()print('影响因子分析矩阵：')print(all_data.corr())model = LinearRegression()model.fit(x, y)a = model.intercept_  # 截距b = model.coef_  # 回归系数return a,b

第四个函数主要实现遥感数据的校正，对比在站点处，校正前和校正后的误差。

#------------------4、校正2013年遥感降水数据 by--------------------
'''
说明：利用建立的校正函数计算2013年遥感卫星降水数据，比较和2013年站点数据的差别
输入：站点降水数据S_prec2013，遥感降水数据T_prec2013
输出：误差（相关系数R2,平均绝对误差MAE）;绘制降水数据对比图
'''
def S_val(ST_prec2013, a, b):ST_prec2013 = np.array(ST_prec2013, dtype=float)s_prec =  ST_prec2013[:,2]                #实际降水即站点降水T_prec =  ST_prec2013[:,3]                #校正前卫星降水x =  ST_prec2013.T[np.arange(3,7)]x=x.Tb=np.array(b)a = a * np.ones(len(ST_prec2013))a=a.T# print(x.shape,b.shape,a.shape)T_revise_2013 = a + np.dot(x,b)                 #校正后卫星降水#误差计算MAE_s = sum(s_prec) / (len(ST_prec2013))  #站点均值MAE_t = sum(T_prec) / (len(T_prec))        #校正前遥感均值MAE_t_revise = sum(T_revise_2013) / (len(T_revise_2013))#校正后遥感均值def val(MAE_s0,MAE_t0,T_revise0):MAE=0R1=0R2 = 0R3 = 0R4=0for i in range(len(ST_prec2013)):MAE = abs(s_prec[i] - T_revise0[i]) + MAER1 = (s_prec[i]-MAE_s0)*(T_revise0[i] - MAE_t0)+R1R2 = (s_prec[i]-MAE_s0)*(s_prec[i]-MAE_s0)+ R2R3 = (T_revise0[i] - MAE_t0) * (T_revise0[i] - MAE_t0) + R3R4 =(s_prec[i] - T_revise0[i])*(s_prec[i] - T_revise0[i])+R4R =R1/(pow(R1*R2,0.5))       #相关系数MAE=MAE/len(ST_prec2013)     #均值RMSE = pow(R4/len(ST_prec2013),0.5)                           #RMSEreturn R,MAE,RMSER1, MAE1,RMSE1 = val(MAE_s,MAE_t,T_prec)R2, MAE2,RMSE2 = val(MAE_s, MAE_t_revise,T_revise_2013)num=np.linspace(0,len(s_prec)-1,num=len(s_prec))ln1,=plt.plot(num, s_prec)ln2, = plt.plot(num, T_prec)ln3,=plt.plot(num, T_revise_2013)plt.title('校准前后降水量数据对比图')  # 折线图标题plt.rcParams['font.sans-serif'] = ['SimHei']  # 显示汉字plt.xlabel('序号')  # x轴标题plt.ylabel('降水量(mm)')  # y轴标题plt.legend(handles=[ln1,ln2,ln3],labels=['s_prec_2013','T_prec_2013','T_revise_2013'],loc='upper right')# plt.show()plt.savefig("02校准前后降水量数据对比图.jpg", dpi=100)plt.clf()return MAE1,R1,RMSE1,MAE2,R2,RMSE2

第五个函数主要实现对空间上遥感数据的整个校正，因为是一定的区域，则需要进行研究区的裁剪和最终结果保存为tiff。该步骤与函数2具有一定的重合，其实可以尝试合并一部分内容，简化函数。

# ------------------5、校正后空间遥感降水数据可视化 by--------------------
'''
说明：利用建立的校正函数计算2013年遥感卫星降水数据
输入：校正前遥感降水数据T_prec2013
输出：校正后的遥感降水数据T_prec2013_corr,最终结果需保存为tiff格式
'''def T_val(DEM_path, T_prec2013_path,a,b):dataset1 = gdal.Open(DEM_path + '\dem_china.tif')dataset2 = gdal.Open(T_prec2013_path + '3B43_2013.TIF')img_geotrans1 = dataset1.GetGeoTransform()img_geotrans2 = dataset2.GetGeoTransform()img_proj2 = dataset2.GetProjection()samples1 = dataset1.RasterXSizelines1 = dataset1.RasterYSizesamples2 = dataset2.RasterXSizelines2 = dataset2.RasterYSizeim_data1 = dataset1.ReadAsArray(0, 0, samples1, lines1)im_data2 = dataset2.ReadAsArray(0, 0, samples2, lines2)del dataset1del dataset2#裁剪图像def clip(pnt_coordinates,img_geotrans):im_loc = [[], []]for i in range(len(pnt_coordinates[0])):sample = sp.Symbol('sample')line = sp.Symbol('line')a = -pnt_coordinates[0][i] + img_geotrans[0] + sample*img_geotrans[1] + line*img_geotrans[2]b = -pnt_coordinates[1][i] + img_geotrans[3] + sample*img_geotrans[4] + line*img_geotrans[5]answer = sp.solve([a,b],[sample,line])#解二元一次方程组line = int(np.floor(answer[line]))sample = int(np.floor(answer[sample]))im_loc[0] += [line]im_loc[1] += [sample]return im_locpnt_coordinates = [[111.25, 121.25], [43, 36]] #左上，右下经纬度im_loc1 = clip(pnt_coordinates, img_geotrans1)im_loc2=clip(pnt_coordinates,img_geotrans2)# print(im_loc1,im_loc2)start_line1 = im_loc1[0][0]end_line1 = im_loc1[0][1]start_sample1 = im_loc1[1][0]end_sample1 = im_loc1[1][1]start_line2 = im_loc2[0][0]end_line2 = im_loc2[0][1]start_sample2 = im_loc2[1][0]end_sample2 = im_loc2[1][1]data1 = im_data1[start_line1:end_line1, start_sample1:end_sample1]data1[data1==32767]=0  #将取值区域的dem去除data2 = im_data2[start_line2:end_line2, start_sample2:end_sample2]#根据校正函数进行校正data2_revise=np.zeros(data2.shape)lon=np.linspace(111.25,121.25,num=len(data2[0]))lat=np.linspace(43,36,num=len(data2))for i in range(len(data2)):for j in range(len(data2[0])):if data1[i,j]==0:data2_revise[i, j]=0else:data2_revise[i,j]=a+b[0]*data2[i,j]+b[3]*data1[i,j]+b[1]*lon[j]+b[2]*lat[i]im_data3=np.zeros([lines2,samples2])im_data3[start_line2:end_line2,start_sample2:end_sample2]=data2_revise #替换计算矩阵块#创建GEOTIFF文件im_width = samples2im_height = lines2im_bands = 1datatype = gdal.GDT_Float32driver = gdal.GetDriverByName("GTiff")dataset = driver.Create("3B43_2013_revise.TIF", im_width, im_height, im_bands, datatype)dataset.SetGeoTransform(img_geotrans2)  # 写入仿射变换参数dataset.SetProjection(img_proj2)  # 写入投影dataset.GetRasterBand(1).WriteArray(im_data3)  # 写入数组数据del  dataset

从箱型图和影响因子分析矩阵中可以看出，站点降水S_prec和遥感提取的降水T_prec、经度lon、纬度lat、所在高程dem均有线性关系，而S_prec和T_prec线性关系属于强相关，符合常识。

（2）校正函数

获得的多元回归函数为：
y= -453.51 + 0.8996 * T_prec + 7.637 * lon-10.9 * lat + 0.04 * dem
其中，y为校正后降水量（mm），T_prec为遥感数据TRMM3B43获得的降水量（mm），lon为经度（度），lat为纬度（度），dem为高程（m）。

（3）误差分析
选取的误差分析评价指标如下。其中，平均绝对误差和均方根误差表示TRMM 3B43数据与实测数据之间的偏差，数值越小越好，相关系数用来衡量 TRMM3B43数据和实测数据之间的相关性，越接近于1表示相关性越好。
平均绝对误差：

均方根误差：

相关系数：


图3 2013遥感数据校正前（左）和校正后（右）对比图

计算2013年遥感数据校正前和校正后，效果如上图，其中，校正后数据，因DEM数据缺失（沿海区域），故右图有部分缺失（黑色圈处）。校正后的影像已保存为tiff格式，可利用ArcGIS直接打开查看效果。

图4 2013年站点位置校正前和校正后对比图

计算2013年站点位置校正前和校正后TRMM3B43数据和站点数据（真值）的误差，平均绝对误差降低，均方根误差降低，相关系数基本保持不变，说明该校正函数校正效果较好。通过与文献（石玉立，2015）结果对比，该方法基本满足校正需求。同时，因训练数据以及验证数据均较少，其实不具有较好的代表性，本方法需要更多的数据验证。

因代码为多位同学共同完成，部分代码习惯未做统一，部分代码存在重复问题，另外，路径问题，本来可以写成相对路径，偷懒就没改，希望后续可以有人补充简化。

部分代码参考：
https://blog.csdn.net/qq_37970770/article/details/103097974
思路参考文献：
石玉立,宋蕾.1998-2012年青藏高原TRMM 3B43降水数据的校准[J].干旱区地理,2015,38(05):900-911.