python根据卫星六根数和两行根数计算卫星经纬度

转载

mob64ca140e76c8 2024-08-02 10:03:03

文章标签 python 开发语言数据预处理换出 文章分类 Python 后端开发

作业要求：

目的：学习用地面站点实测降雨数据校正卫星遥感降雨数据。
原始数据：
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年时间序列数据。

详细过程

文章目录

1、总体思路
2、代码详解
3、结果分析
后记

1、总体思路

选取2009-2013年的年降水数据，对京津冀及周边地区的大概范围（经度：111.25-121.25，纬度：36-43）进行影像降水数据的校正。前提假设站点获得的降水数据为真值。另外，增加了全国dem数据，提高校正精度。总体思路如下：

（1）读取2009-2013年，研究区内各站点的年降水量、经度、纬度、DEM、站点所处位置对应的遥感数据TRMM3B43观测得到的降水量。

（2）利用2009-2012年的数据，分析站点的年降水量和经度、纬度、DEM、遥感数据TRMM3B43观测得到的降水量之间的相关性，同时，建立多元回归函数。

（3）利用2013年的遥感数据TRMM3B43观测得到的降水量，经度、纬度、DEM计算校正后的降水量，该降水量与2013年站点降水数据比较，进行误差分析，评估校正效果。

python根据卫星六根数和两行根数计算卫星经纬度_数据

图1 总体技术路线

2、代码详解

主函数，将整个过程分为5步完成，包含五个函数来实现。

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 LinearRegression

if __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_line
        stime = stime.strip()
        if stime in years:
            try:
                # 做一个判断，剔除降水值为负的值
                i_prec = int(sprec.strip())
                if i_prec > 0:
                    return stime.strip(), sname.strip(), i_prec
            except Exception:
                return "", "", 0
        return "", "", 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.1
                s_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.1
                s_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.RasterXSize
    lines = dataset1.RasterYSize
    # bands = dataset1.RasterCount
    img_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 dataset1
    del dataset2
    del dataset3
    del dataset4
    del 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.T

    all_data =  S_prec200912.T[np.arange(2,7)]
    all_data=all_data.T
    all_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.T
    b=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=0
        R1=0
        R2 = 0
        R3 = 0
        R4=0
        for i in range(len(ST_prec2013)):
            MAE = abs(s_prec[i] - T_revise0[i]) + MAE
            R1 = (s_prec[i]-MAE_s0)*(T_revise0[i] - MAE_t0)+R1
            R2 = (s_prec[i]-MAE_s0)*(s_prec[i]-MAE_s0)+ R2
            R3 = (T_revise0[i] - MAE_t0) * (T_revise0[i] - MAE_t0) + R3
            R4 =(s_prec[i] - T_revise0[i])*(s_prec[i] - T_revise0[i])+R4
        R =R1/(pow(R1*R2,0.5))       #相关系数
        MAE=MAE/len(ST_prec2013)     #均值
        RMSE = pow(R4/len(ST_prec2013),0.5)                           #RMSE
        return R,MAE,RMSE
    R1, 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.RasterXSize
    lines1 = dataset1.RasterYSize
    samples2 = dataset2.RasterXSize
    lines2 = dataset2.RasterYSize
    im_data1 = dataset1.ReadAsArray(0, 0, samples1, lines1)
    im_data2 = dataset2.ReadAsArray(0, 0, samples2, lines2)
    del dataset1
    del 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_loc

    pnt_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]=0
                else:
                    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 = samples2
    im_height = lines2
    im_bands = 1
    datatype = gdal.GDT_Float32

    driver = 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