python3 用地面站点实测降雨数据校正TRMM/3B43遥感降雨数据

目的：学习用地面站点实测降雨数据校正卫星遥感降雨数据。
学习将多景遥感图像绘制成在一张图上，且使用同一个colorbar
原始数据：
2007-2010年的TRMM卫星3B43产品（覆盖范围精度-180~ 180、纬度-50~50，0.25x0.25度分辨率，像素值单位：mm/hour）；
江西梅川江地区63个气象站的2007-2010年的月降雨数据，以及站点位置SHAPE文件；
数据处理要求：
实验区2（江西梅川江）
对覆盖范围（经度：115-117,纬度：25.5-27.5）内的卫星年降雨数据，用地面数据进行校正，并存为TIF文件。
校正的卫星年降雨数据可以是一年的，如：2010年，也可以是2007-2010年时间序列数据。
两个实验区，选择一组数据处理即可
建议采用python和本课程学习的各种软件包完成，但是本着“抓着老鼠就是好猫”原则，可不局限于本门课程讲授的内容。

引言

在遥感领域处理方法即是将站点附近3*3像元内的平均值作为遥感图像站点值与地面站点值进行比较，或以站点所在位置为圆心，10km半径（或其他半径）内的像元平均值作为遥感图像站点值与地面站点值进行比较，选择1/3的实验数据用来建模，2/3的实验数据用来验证模型。

在GIS(Geographic Information Science,GIS)领域中，由于尺度问题造成的误差是通过另一种尺度匹配方法减小的，通常对站点数据shape建立泰森多边形，将点尺度转换为栅格尺度，每个栅格值基于不同因素加权平均得到，然后将转换后的栅格像元与遥感图像像元值进行比较建模研究，选择1/3的实验数据用来建模，2/3的实验数据用来验证模型。

本实验采用前者，但考虑到TRMM像元空间分辨率为2.5km，63个站点之间的距离在TRMM遥感图像上非常接近，无法取3*3像元或圆心半径法进行研究。本实验直接取2010年地面站点值与对应位置处遥感图像值进行建模，由于梅江川地区63个雨量站位置相近，海拔、经纬度等因素对其模型影响可忽略不计，因此只考虑数值误差之间的差异，选择一元线性回归方法进行建模。

建立模型

点击链接，查看上一篇文章。
用地面雨量计数据校正TRMM/3B43数据之模型建立

模型应用

使用上述线性模型进行校正，并将校正前后的梅川江地区遥感图像绘制成在一张图上，且使用同一个colorbar。
前篇代码
改进后的代码，如下：

# -*- coding: utf-8 -*-
"""
Created on Sat Nov 16 09:43:08 2019

@author: Administrator
"""

from osgeo import gdal
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

#-------------------------读取2007、2008、2009年TRMM降水遥感图像-----------------------
path = 'E:\\课程PPT\\Python空间数据处理\\python空间数据处理-期中大作业02\\卫星降雨数据TRMM3B43_分辨率0.25度\\'
#有时候路径必须双反斜杠
dataset = gdal.Open(path + '3B43_2007.TIF')
dataset2 = gdal.Open(path + '3B43_2008.TIF')
dataset3 = gdal.Open(path + '3B43_2009.TIF')
samples = dataset.RasterXSize
lines = dataset.RasterYSize
bands = dataset.RasterCount
img_geotrans = dataset.GetGeoTransform()
#Out[10]: (-180.0, 0.25, 0.0, 50.0, 0.0, -0.25)
img_proj = dataset.GetProjection()
#Out[12]: 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84"...]]]
#可以看到此遥感图像未进行投影定义，因此还不是投影坐标，只是地理坐标
#在求解这些实测站点位于图像中哪个行列位置时不用将其转为投影坐标
im_data1 = dataset.ReadAsArray(0,0,samples,lines)
im_data2 = dataset2.ReadAsArray(0,0,samples,lines)
im_data3 = dataset3.ReadAsArray(0,0,samples,lines)
del dataset
del dataset2
del dataset3
#--------------------------------校正数据---------------------------------------
cor_data1 = 0.797 * im_data1 + 352.279
cor_data2 = 0.797 * im_data2 + 352.279
cor_data3 = 0.797 * im_data3 + 352.279
#-------------------------------找到研究区范围的行列号-------------------------------------
#-----------------------转换左上、右下经纬度坐标到对应遥感图像上的行列号----------
im_loc = [[],[]]
#第一个列表放行号，存纬度换出来的值，后者为列号，存经度换出来的值
# variable += [value] 等同于 variable.append(value)
#左上角 115 27.5
#右下角 117 25.5
pnt_coordinates = [[115,117],[27.5,25.5]]
for i in range(len(pnt_coordinates[0])):
    sample = sp.Symbol('sample')
    line = sp.Symbol('line')
    #pnt_coordinates[0][i] = img_geotrans[0] + sample*img_geotrans[1] + line*img_geotrans[2]
    #pnt_coordinates[1][i] = img_geotrans[3] + sample*img_geotrans[4] + line*img_geotrans[5]
    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])#解二元一次方程组
    #Out[15]: {x: 0, y: 1}
    line = int(np.floor(answer[line]))
    sample = int(np.floor(answer[sample]))
    im_loc[0] += [line]
    im_loc[1] += [sample]
#------------------------------起始终止行列号------------------------------------
start_line = im_loc[0][0]
end_line = im_loc[0][1]
start_sample = im_loc[1][0]
end_sample = im_loc[1][1]
#将所有裁剪的数据存在一个列表中
datasets = []
for j in range(3):
    datasets += [eval('im_data'+str(j+1))[start_line:end_line,start_sample:end_sample]]
for i in range(3):
    datasets += [eval('cor_data'+str(i+1))[start_line:end_line,start_sample:end_sample]]
#----------------------------用于生产colorbar的vmin和vmax-----------------------
vmin = np.min(datasets)
vmax = np.max(datasets)
#------------------------------------绘图---------------------------------------
fig,axs = plt.subplots(nrows = 2, ncols = 3,figsize = (14,8))
#figsize = (width，hight) #按照出图比例来
extent = (0,1,0,1)
#将x周和y轴都归一化。方便指定label的位置(0,1)相对比例，否则指定label位置时不好指定
norm = Normalize(vmin = vmin,vmax = vmax)
#Normalize()跟归一化没有任何关系，函数的作用是将颜色映射到vmin-vmax上，即让接下来的颜色表/颜色柱的起始和终止分别取值vmin和vmax
for i in range(2):
    if i == 0:        
        for j in range(3):
            axs[i,j].imshow(datasets[j],extent = extent,norm = norm,cmap = 'jet_r')
            axs[i,j].set_xlabel(str(2007+j)+ ' Original',fontdict = {'family' : 'Times New Roman','fontweight':'bold','size':15})
            axs[i,j].set_xticks(np.linspace(0,1,5))
            font = {'fontweight':'bold'}
            axs[i,j].set_xticklabels(('115.0E','115.5E','116.0E','116.5E','117.0E'),fontdict = font)
            axs[i,j].set_yticks(np.linspace(0,1,5))
            axs[i,j].set_yticklabels(('25.5 N','26.0 N','26.5 N','27.0 N','27.5 N'),fontdict = font)
    else:
        for j in range(3):
            axs[i,j].imshow(datasets[j+3],extent = extent,norm = norm,cmap = 'jet_r')
            axs[i,j].set_xlabel(str(2007+j)+ ' Corrected',fontdict = {'family' : 'Times New Roman','fontweight':'bold','size':15})
            axs[i,j].set_xticks(np.linspace(0,1,5))
            font = {'fontweight':'bold'}
            axs[i,j].set_xticklabels(('115.0E','115.5E','116.0E','116.5E','117.0E'),fontdict = font)
            axs[i,j].set_yticks(np.linspace(0,1,5))
            axs[i,j].set_yticklabels(('25.5 N','26.0 N','26.5 N','27.0 N','27.5 N'),fontdict = font)
            if j == 2:
                sc = axs[i,j].imshow(datasets[j+3],extent = extent,norm = norm,cmap = 'jet_r')
#一行三个子图的总宽度 为 全部宽度的 0.9；剩下的0.1用来放置colorbar
fig.subplots_adjust(right=0.9)
#colorbar 左 下 宽 高 
l = 0.92
b = 0.53
w = 0.015
h = 0.35
#对应 l,b,w,h；设置colorbar位置；
rect = [l,b,w,h] 
cbar_ax = fig.add_axes(rect) 
plt.colorbar(sc, cax=cbar_ax)
plt.savefig('correct_pictures.png',dpi = 300,bbox_inches='tight')
plt.show()

模型精度评价

为评价线性模型校正精度，分别计算校正之前和校正之后的TRMM/3B43遥感反演数据与雨量计数据之间的MAE(mean absolute error,MAE)、RMSE(root mean square error,RMSE)和 ，计算63个雨量计数据与校正前后的上述统计指标。
注意：接下来的代码需在上面的基础上运行

# -*- coding: utf-8 -*-
"""
Created on Tue Nov 12 19:31:07 2019

@author: BAI Depei
"""

from osgeo import ogr
from osgeo import gdal
import numpy as np
import sympy as sp

#------------------read the station_coordinates--------------------
fn = 'E:\课程PPT\Python空间数据处理\python空间数据处理-期中大作业02\梅川江地面数据\meichuan_prec_station.shp'
ds = ogr.Open(fn,False)
layer = ds.GetLayer(0)
spatialref = layer.GetSpatialRef()
lydefn = layer.GetLayerDefn()  #图层定义信息
geomtype = lydefn.GetGeomType()
fieldlist = [] 
for i in range(lydefn.GetFieldCount()):
    fddefn = lydefn.GetFieldDefn(i)
    fddict = {'name':fddefn.GetName(),'type':fddefn.GetType(),\
              'width':fddefn.GetWidth(),'decimal':fddefn.GetPrecision()}
    fieldlist += [fddict]
    #获取属性表每列列名及其属性：类型、宽度、精度
    geomlist,reclist = [],[]
    #geomlist存放所有图形的坐标点(字符串中包含)，reclist存放所有的图形名
    feature = layer.GetNextFeature()
    #首次获取feature
    while feature is not None:
        geom = feature.GetGeometryRef()
        geomlist += [geom.ExportToWkt()]
        rec = {}
        for fd in fieldlist:
            rec[fd['name']] = feature.GetField(fd['name'])
        reclist += [rec]
        feature = layer.GetNextFeature()
ds.Destroy()
#关闭数据
#-------------------------将坐标写入列表中---------------------------------
pnt_coordinates = [[],[]]
for i in range(len(geomlist)):
    a = geomlist[i].split(' ')
    b = a[1].split('(')
    c = a[2].split(')')
    longitude = float(b[1])
    latitude = float(c[0])
    pnt_coordinates[0] += [longitude]
    pnt_coordinates[1] += [latitude]
#reclist则不必写入其他变量，因其调用简单

#-------------------------读取2010年TRMM降水遥感图像-----------------------
path = 'E:\\课程PPT\\Python空间数据处理\\python空间数据处理-期中大作业02\\卫星降雨数据TRMM3B43_分辨率0.25度\\'
#有时候路径必须双反斜杠
dataset = gdal.Open(path + '3B43_2010.TIF')
samples = dataset.RasterXSize
lines = dataset.RasterYSize
bands = dataset.RasterCount
img_geotrans = dataset.GetGeoTransform()
#Out[10]: (-180.0, 0.25, 0.0, 50.0, 0.0, -0.25)
img_proj = dataset.GetProjection()
#Out[12]: 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84"...]]]
#可以看到此遥感图像未进行投影定义，因此还不是投影坐标，只是地理坐标
#在求解这些实测站点位于图像中哪个行列位置时不用将其转为投影坐标
im_data = dataset.ReadAsArray(0,0,samples,lines)

del dataset
#-----------------------转换站点经纬度坐标到对应遥感图像上的行列号----------
im_loc = [[],[]]
#第一个列表放行号，存纬度换出来的值，后者为列号，存经度换出来的值
# variable += [value] 等同于 variable.append(value)
for i in range(len(pnt_coordinates[0])):
    sample = sp.Symbol('sample')
    line = sp.Symbol('line')
    #pnt_coordinates[0][i] = img_geotrans[0] + sample*img_geotrans[1] + line*img_geotrans[2]
    #pnt_coordinates[1][i] = img_geotrans[3] + sample*img_geotrans[4] + line*img_geotrans[5]
    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])#解二元一次方程组
    #Out[15]: {x: 0, y: 1}
    line = int(np.floor(answer[line]))
    sample = int(np.floor(answer[sample]))
    im_loc[0] += [line]
    im_loc[1] += [sample]
#----------------------read 2007_2009_TRMM_stationpnt_value-------------------------    
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] += [cor_data1[im_loc[0][i],im_loc[1][i]]]
    data_TRMM_stationpnt[4] += [cor_data2[im_loc[0][i],im_loc[1][i]]]
    data_TRMM_stationpnt[5] += [cor_data3[im_loc[0][i],im_loc[1][i]]]
#----------------------read 2010_rain_gauge_value------------------------------
path2 = 'E:\\课程PPT\\Python空间数据处理\\python空间数据处理-期中大作业02\\梅川江地面数据\\'
file_handle = open(path2 + 'meichuan_precipitation_2007-2010.txt')
line_value = file_handle.readlines()
file_handle.close()

#data_TRMM_stationpnt是根据geomlist次序依次遍历获得的，
#而geomlist次序与reclist次序一致，因此用reclist中台站名和2010作为关键字去检索
data_rain_gauge_2007_2009 = [[],[],[]]
#站点数据
for i in range(1,len(line_value)):
    temp = line_value[i][:-2].split(',')
    #Out[3]: '200704,DongShao,184.200000\n'
    if temp[0] == '200701':
        data = []
        for j in range(12):
            temp  = line_value[i+j][:-2].split(',')
            data += [float(temp[2])]
        total_value = sum(data)
        data_rain_gauge_2007_2009[0] +=  [total_value]
        i += 12
        continue
    if temp[0] == '200801':
        data = []
        for j in range(12):
            temp  = line_value[i+j][:-2].split(',')
            data += [float(temp[2])]
        total_value = sum(data)
        data_rain_gauge_2007_2009[1] +=  [total_value]
        i += 12
        continue
    if temp[0] == '200901':
        data = []
        for j in range(12):
            temp  = line_value[i+j][:-2].split(',')
            data += [float(temp[2])]
        total_value = sum(data)
        data_rain_gauge_2007_2009[2] +=  [total_value]
        i += 12
        continue
#------------------------calculate statistics index--------------------------
mean_TRMM_2007 = np.mean(data_TRMM_stationpnt[0])
mean_TRMM_2008 = np.mean(data_TRMM_stationpnt[1])
mean_TRMM_2009 = np.mean(data_TRMM_stationpnt[2])
MAE = []
for i in range(3):
    MAE += [sum(abs(np.array(data_TRMM_stationpnt[i])-np.array(data_rain_gauge_2007_2009[i])))/63]
for i in range(3):
    MAE += [sum(abs(np.array(data_TRMM_stationpnt[i+3])-np.array(data_rain_gauge_2007_2009[i])))/63]
RMSE = []
for i in range(3):
    RMSE += [np.sqrt((sum(abs(np.array(data_TRMM_stationpnt[i])-np.array(data_rain_gauge_2007_2009[i])))**2)/63)]
for i in range(3):
    RMSE += [np.sqrt((sum(abs(np.array(data_TRMM_stationpnt[i+3])-np.array(data_rain_gauge_2007_2009[i])))**2)/63)]
R = []
for i in range(6):
    if i <3:       
        mean_TRMM = np.mean(np.array(data_TRMM_stationpnt[i]))
        mean_rain = np.mean(np.array(data_rain_gauge_2007_2009[i]))
        a = np.array(data_TRMM_stationpnt[i]) - mean_TRMM
        b = np.array(data_rain_gauge_2007_2009[i]) - mean_rain
        c = (sum(a*b))/np.sqrt(sum(a**2)*sum(b**2))
        R.append(c)
        i += 1 
    else:
        mean_TRMM = np.mean(np.array(data_TRMM_stationpnt[i]))
        mean_rain = np.mean(np.array(data_rain_gauge_2007_2009[i-3]))
        a = np.array(data_TRMM_stationpnt[i]) - mean_TRMM
        b = np.array(data_rain_gauge_2007_2009[i-3]) - mean_rain
        c = (sum(a*b))/np.sqrt(sum(a**2)*sum(b**2))
        R += [c]

结论

结论就是该模型效果不好，此处不予展示。
通过该模型校正之后，TRMM数据与雨量计实测数据之间的MAE明显变大，RMSE相比校正前明显变大，模型校正前后R变化不大，因此可知利用2010年建立的线性回归模型虽然本身是显著的，但利用该模型校正其他年份的TRMM数据效果不理想，考虑到各个雨量站之间的地理环境差异不大，因此建立与地理环境因素有关的模型是不成立的，而纯粹的一元线性回归数学模型又会使得偏差增大，因此下一步可以考虑尝试其他数学模型如残差模型等来进行校正。

---

好好生活，珍惜爱你的和你爱的人!

版权归作者 小白是哪个小白_ 所有，转载、引用请注明链接出处。