python3 用地面站点实测降雨数据校正TRMM/3B43遥感降雨数据 |
您所在的位置：网站首页 › trmm数据介绍 › python3 用地面站点实测降雨数据校正TRMM/3B43遥感降雨数据 |
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 |
【本文地址】
今日新闻 |
推荐新闻 |