一、不同差值效果对比

二、制图代码

2.1用到的模块

2.1.1遇到的问题

2.2经纬度转shp坐标点

2.2.1遇到的问题

2.3IDW(反距离权重)

2.3.1遇到的问题

2.4利用Basemap渲染

2.4.1遇到的问题

三、结果展示

3.1arcpy效果

3.2Basemap模块效果

一、不同差值效果对比

因为arcpy不能部署在linux因此采用python第三方模块，左侧为反距离权重差值方法，右侧为最近邻，等方法，可以直观的看出空间插值，idw效果较好。

二、制图代码

2.1用到的模块

因为最终程序需要部署到linux系统服务器上，arcgis等工具只能用于windows下，因此采用python第三方包。

import os
import conda
conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
proj_lib = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib
#以上代码是linux中运行可能产生报错，需要加上
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import cx_Oracle
import shapefile as sf
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import datetime
from matplotlib.font_manager import FontProperties
#windows下需要定义字体中文
font = FontProperties(fname='C:/Windows/Fonts/simsun.ttc', size=13)import gdal
import osgeo.ogr as ogr
import osgeo.osr as osr
#linux下没有gpu的做法。
plt.switch_backend('agg')

2.1.1遇到的问题

头部导入模块代码，其中除了import，很多内容是为了防止报错。

错误一，可能是因为环境问题，安装包后会报类似No such file or directory的错误，因此需要加上

conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
proj_lib = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib

错误二，matplotlib不能支持中文，试过了很多方法，指定字体文件，window和linux下通用，字体可以自己选择，只需要使用的时候，指定font，代码如下。

font = FontProperties(fname='C:/Windows/Fonts/simsun.ttc', size=13)

错误三，linux下plt模块展示可能报错，需要指定用那种方式，代码如下，参数可以自行搜素。

plt.switch_backend('agg')

2.2经纬度转shp坐标点

gdal提供的差值函数，参数为shp点数据和最终的输出结果，因此这里将读取到的经纬度数据转换为shp，当然也可选择自己写反距离权重插值的算法，不将数据转换shp点，进行插值。

def XY_Point(shp_path,data_list):""":param shp_path: 输出文件路径:param data_list: 包含x，y，data的数据集[[112.503, 34.844, 22.0], [112.09, 34.43, 21.0]]:return:"""#定义驱动driver = ogr.GetDriverByName("ESRI Shapefile")#创建一个shp文件data_source = driver.CreateDataSource(shp_path)srs = osr.SpatialReference()srs.ImportFromEPSG(4326) #这是WGS84,想用什么自己去搜下对应的编码就行了#创建点文件layer = data_source.CreateLayer("Point", srs, ogr.wkbPoint)#定义字段的内容，添加字段field_name = ogr.FieldDefn("data", ogr.OFTString)field_name.SetWidth(14)layer.CreateField(field_name)feature = ogr.Feature(layer.GetLayerDefn())for xx in data_list:x = xx[0]y = xx[1]data = xx[2]#写入数据feature.SetField("data", "{0}".format(str(data)))#生成点的固定格式wkt = "POINT({0} {1})".format(x, y)point = ogr.CreateGeometryFromWkt(wkt)feature.SetGeometry(point)layer.CreateFeature(feature)feature = Nonedata_source = None

2.2.1遇到的问题

问题一，坐标问题，坐标简单的可以分为地理坐标系和投影坐标系，经纬度为地理坐标系，WGS84，编号为4326，坐标系一定要注意，因为涉及到不同模块，默认坐标参数可能不一致，使用中应该调节为一致。

srs.ImportFromEPSG(4326)

问题二，生成点数据的固定格式，矢量数据分为点、线、面、生成矢量的时候，每一种数据都有自己对应的格式，具体可以参考此链接https://blog.csdn.net/weixin_42464154/article/details/104112786。

问题三，怎么将多个坐标点生成在同一个shp文件中，因为网上大多教程都是单一点生成，循环遍历，然后将最后两行代码放到循环外面。

feature = None
data_source = None

2.3IDW(反距离权重)

gdal自带模块，具体度娘有详细的解释。

def idw(output_file, point_station_file):"""idw空间插值:param output_file:插值结果:param point_station_file: 矢量站点数据:return:"""# 代码调用outputBounds定义范围opts = gdal.GridOptions(algorithm="invdistnn:power=2.0:smothing=0.0:radius=1.0:max_points=12:min_points=0:nodata=-10",format="GTiff", outputType=gdal.GDT_Float32, zfield="data",outputBounds=[110,31,117,37])gdal.Grid(destName=output_file, srcDS=point_station_file, options=opts)

2.3.1遇到的问题

问题一，怎么确定插值的范围，如果范围不能确定，Basemap中的范围虽然可以根据差值的栅格动态调整出图范围，解决坐标偏移问题，但是出图可能裁剪掉感兴趣区域，如下图所示，因此这里的范围要和Basemap裁剪的范围保持一致，调节参数中的范围即可。

outputBounds=[110,31,117,37]

2.4利用Basemap渲染

basemap这个模块确实好用，但是接触的较少，具体用法可以搜索具体开发手册一类的文件。

def chutu(name,tif):""":param name: 判断输入的数据是温度还是降水:param tif: 经过gdal差值后的tif文件路径:return:"""cont=""ds = gdal.Open(tif)grid_z = ds.ReadAsArray()adfGeoTransform = ds.GetGeoTransform()shapexy = grid_z.shape#根据分辨率和左上角坐标计算范围xmx = adfGeoTransform[0] + adfGeoTransform[1] * shapexy[0]ymi = adfGeoTransform[3] + adfGeoTransform[-1] * shapexy[1]#plt.rcParams['font.sans-serif']=['SimHei']# 用来正常显示负号#plt.rcParams['axes.unicode_minus'] = Falsemap = Basemap(projection = 'cyl',llcrnrlon=adfGeoTransform[0], llcrnrlat=adfGeoTransform[3], urcrnrlon=xmx, urcrnrlat=ymi,resolution = 'h')#cyl,projection = 'tmerc'#llcrnrlon=110, llcrnrlat=31, urcrnrlon=117, urcrnrlat=37,map.readshapefile(hn,'comarques',drawbounds = True)#制作河南的图形for info, shape in zip(map.comarques_info, map.comarques):lon, lat = zip(*shape)map.plot(lon, lat, marker = None, color= 'k',lineStyle='-',linewidth=0.5)#--虚线，-.省界map.plot(112.356, 33.476, marker = "o", color= 'r',markersize=6)#南召位置#南召标注plt.annotate('南召县', xy=(112.356, 33.476),  xycoords='data',xytext=(5, 0), textcoords='offset points',arrowprops=None,fontproperties=font,size=10)#grid_x, grid_y = np.mgrid[110:120,31:40]# 插值方法：'nearest', 'linear', 'cubic'np_data = np.array(DATA_list)np_w = np.array(W_list)#坐标点集合#grid_z = griddata(np_w, np_data, (grid_x, grid_y), method='nearest')x1 = np.linspace(map.llcrnrlon, map.urcrnrlon, shapexy[1])y1 = np.linspace(map.llcrnrlat, map.urcrnrlat, shapexy[0])grid_x, grid_y = np.meshgrid(x1, y1)if name=="TEM":cont = '南召县逐小时温度\n %s年%s月%s日%s时BJT' % (year, month, day, hour)plt.title(cont, fontproperties=font, size=13)levels = np.linspace(-6, 38, 50)cbar_levels = np.linspace(-6, 38, 20)# cf = map.scatter(np_x, np_y, np_data, c=np_data, cmap='jet', alpha=0.75)cf = map.contourf(grid_x, grid_y, grid_z, levels=levels, cmap=plt.cm.jet)cbar = map.colorbar(cf, location='right', format='%d℃', size=0.1,ticks=cbar_levels)cbar.ax.tick_params(labelsize=8)  # 图例宽度if name=="RAIN":cont = '南召县日累计降水分布\n %s年%s月%s日BJT' % (year, month, day)plt.title(cont,fontproperties=font,size=13)levels = np.linspace(0, 50, 50)cbar_levels = [1,10,25,50,100,200]#cf = map.scatter(np_x, np_y, np_data, c=np_data, cmap='jet', alpha=0.75)# Colors是一些自选颜色列表Colors = ('#8FF041', '#38A800', '#00A9E6', '#005CE6','#E600A9')#'#CEFFCE', '#28FF28', '#007500', '#FFFF93', '#8C8C00', '#FFB5B5',# '#FF0000', '#CE0000', '#750000')cf = map.contourf(grid_x, grid_y, grid_z, levels=cbar_levels,colors=Colors) #norm=mpl.colors.Normalize(vmin=0, vmax=100))#cmap= plt.cm.cool)cbar = map.colorbar(cf, location='right', format='%dmm', size=0.1,ticks=cbar_levels)cbar.ax.tick_params(labelsize=8)#图例宽度#删除边界外面的数据sjz = sf.Reader(hn)for shape_rec in sjz.shapeRecords():vertices = []codes = []pts = shape_rec.shape.pointsprt = list(shape_rec.shape.parts) + [len(pts)]for i in range(len(prt) - 1):for j in range(prt[i], prt[i+1]):vertices.append(map(pts[j][0], pts[j][1]))codes += [Path.MOVETO]codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)codes += [Path.CLOSEPOLY]clip = Path(vertices, codes)clip = PathPatch(clip, transform=ax.transData)for contour in cf.collections:contour.set_clip_path(clip)plt.savefig('%s_%s.png'%(currenttime,name))

2.4.1遇到的问题

问题一，如何自定义数值以及对应颜色。比如数值1-5位黄色，数值5-10为红色，主要代码如下，修改colors和levels的参数为自己想要数值以及对应的颜色。

cbar_levels = [1,10,25,50,100,200]
Colors = ('#8FF041', '#38A800', '#00A9E6', '#005CE6','#E600A9')
cf = map.contourf(grid_x, grid_y, grid_z, levels=cbar_levels,colors=Colors)

问题二，怎么保证行政区范围外的数据为空值，具体利用了shapefile模块,这里裁剪代码是来自互联网copy，用的话复制粘贴就行了，最终效果如下图所示。

sjz = sf.Reader(hn)for shape_rec in sjz.shapeRecords():vertices = []codes = []pts = shape_rec.shape.pointsprt = list(shape_rec.shape.parts) + [len(pts)]for i in range(len(prt) - 1):for j in range(prt[i], prt[i+1]):vertices.append(map(pts[j][0], pts[j][1]))codes += [Path.MOVETO]codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)codes += [Path.CLOSEPOLY]clip = Path(vertices, codes)clip = PathPatch(clip, transform=ax.transData)for contour in cf.collections:contour.set_clip_path(clip)

三、结果展示

3.1arcpy效果

arcig作为GIS专业工具，制图效果一流（哈哈）。

3.2Basemap模块效果

左侧为降水，右侧为温度，虽然和arcgis存在一定差距，但是总体上还算过的去。

GDAL+Basemap+IDW(反距离权重)代替ARCPY，制作温度、降雨分布图相关推荐

【数据可视化应用】Python反距离权重(IDW)插值计算及可视化绘制
本文我们将介绍IDW(反距离加权法(Inverse Distance Weighted)) 插值的Python计算方法及插值结果的可视化绘制过程.主要涉及的知识点如下: IDW简介自定义Python ...
ArcGIS空间插值方法反距离权重法（IDW）的工作原理
反距离权重 (IDW) 插值使用一组采样点的线性权重组合来确定像元值.权重是一种反距离函数.进行插值处理的表面应当是具有局部因变量的表面. 此方法假定所映射的变量因受到与其采样位置间的距离的影响而减小 ...
python反距离权重法_使用Python进行反距离加权(IDW)插值
10月20日改变:这个类Invdisttree组合了反距离权重 scipy.spatial.KDTree. 忘记原来的强力回答; 这是分散数据插值的选择方法. """ i ...
python实现反距离权重插值(IDW)
原文:https://mp.weixin.qq.com/s/2y13vGtqj55Fae-72YXY6g 1 什么叫反距离权重插值? 反距离权重:距离未知点最近的点分配的权重较大,且权重作为距离的函数 ...
反距离权重插值（IDW）
反距离权重插值适用于整个研究区数据均匀分布,不存在聚类的情况,效果最优! 插值步骤步骤一:启动地统计向导,选择确定性方法中的反距离权重法,选择插值的源数据集和数据字段,如有必要,可添加权重字段,完毕 ...
反距离权重加权插值的理解及Python实现
IDW(反距离里加权插值) 假设距离较近的事物要比距离较远的事物更相似. 当为任何未测量的位置预测值时,反距离权重法会采用预测位置周围的测量值与距离预测位置较远的测量值相比,距离预测位置最近的测量值 ...
python反距离权重法_反距离权重法 (Spatial Analyst)—ArcMap | 文档
使用反距离权重法 (IDW) 获得的像元输出值限定在插值时用到的值范围之内.因为反距离权重法是加权平均距离,所以该平均值不可能大于最大输入或小于最小输入.因此,如果尚未对这些极值采样,便无法创建山脊或 ...
gstat | 空间插值（一）——反距离权重插值；使用ggplot2绘制地图
本篇既是空间插值系列的第一篇推文,也是ggplot2工具包系列推文中的一篇.空间插值使用的工具包是gstat,该工具包主要用于地统计分析. library(gstat) 示例数据来自HSAR工具包: ...
反距离加权插值法例题_GMS插值中的反距离权重法（Inverse distance weighted）
反距离权重法是GMS地层插值的默认方法,了解一些关于它的原理会帮助得到更好的插值结果.这次主要介绍Shepard's method方法.反距离权重法基本思路:插值点受附近点的影响最大,而距离较远的点的 ...

GDAL+Basemap+IDW(反距离权重)代替ARCPY，制作温度、降雨分布图

一、不同差值效果对比

二、制图代码

2.1用到的模块

2.1.1遇到的问题

2.2经纬度转shp坐标点

2.2.1遇到的问题

2.3IDW(反距离权重)

2.3.1遇到的问题

2.4利用Basemap渲染

2.4.1遇到的问题

三、结果展示

3.1arcpy效果

3.2Basemap模块效果

GDAL+Basemap+IDW(反距离权重)代替ARCPY，制作温度、降雨分布图相关推荐

最新文章

热门文章