GDAL+Basemap+IDW(反距离权重)代替ARCPY，制作温度、降雨分布图
目录一、结果展示1.1arcpy效果二、制图代码2.1用到的模块2.1.1遇到的问题2.2经纬度转shp坐标点2.2.1遇到的问题2.3IDW(反距离权重)2.3.1遇到的问题2.4利用Basemap渲染2.4.1遇到的问题一、结果展示因为arcpy不能部署在linux因此采用python第三方模块，左侧为反距离权重差值方法，右侧为最近邻，等方法，可以直观的看出空间插值，idw效果较好。1.1ar
·
目录
一、不同差值效果对比
因为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 = None
data_source = None2.2.1遇到的问题
问题一，坐标问题，坐标简单的可以分为地理坐标系和投影坐标系，经纬度为地理坐标系，WGS84，编号为4326，坐标系一定要注意 ，因为涉及到不同模块，默认坐标参数可能不一致，使用中应该调节为一致。
srs.ImportFromEPSG(4326)问题二，生成点数据的固定格式，矢量数据分为点、线、面、生成矢量的时候，每一种数据都有自己对应的格式，具体可以参考此链接https://blog.csdn.net/weixin_42464154/article/details/104112786。
问题三，怎么将多个坐标点生成在同一个shp文件中，因为网上大多教程都是单一点生成，循环遍历，然后将最后两行代码放到循环外面。
feature = None
data_source = None2.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'] = False
map = 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.points
prt = 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.points
prt = 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存在一定差距，但是总体上还算过的去。
更多推荐
4
0- 0
已为社区贡献1条内容
所有评论(0)