哈喽、哈喽大家好！！
最近浏览了不少代码，get到不少新的知识！！！
接下来就直接给大家分享一下，有需要的小伙伴直接打包带走就好了！

文章目录

前言
库函数准备
分段讲解
- 添加比例尺
- 添加比例尺
- 图像裁剪代码
- 栅格数据读取
- 画图函数
完整代码展示
博主说两句

前言

最近用GIS在批量出图，发现一张一张出图真的麻烦（那个累啊！！！）
于是便有了今天这篇文章，初步教大家如何用Python出那种开起来专业一点点的GIS图。

库函数准备

本次使用的库函数也不是很多主要就是以下的这些

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as cor
import cartopy.io.shapereader as sr
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import cartopy.crs as ccrs
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapefile

分段讲解

声明：添加比例尺、添加指北针等并非我原创，也都是搜集于CSDN、Github等等等之类的
但是我或多或少对原本的代码进行了修改，使之更好用一些。

添加比例尺

原始代码中包含了三种样式的图例，样子都还不错。
ax：是我们创建的子图
lon，lat：是我们图例想要放在那个位置的坐标，根据个人喜好来！！！
length：是我们比例的你所输入的比例，比如200等
size：是控制比例尺的高度的（比例尺上三根竖线的高度，一会下面会有展示的）

#-----------函数：添加比例尺--------------
def add_scalebar(ax,lon0,lat0,length,size=0.45):'''ax: 坐标轴lon0: 经度lat0: 纬度length: 长度size: 控制粗细和距离的'''# style 3ax.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=1, label='%d km' % (length))ax.vlines(x = lon0, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)ax.vlines(x = lon0+length/2/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)ax.vlines(x = lon0+length/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)ax.text(lon0+length/111,lat0+size+0.05,'%d' % (length),horizontalalignment = 'center')ax.text(lon0+length/2/111,lat0+size+0.05,'%d' % (length/2),horizontalalignment = 'center')ax.text(lon0,lat0+size+0.05,'0',horizontalalignment = 'center')ax.text(lon0+length/111/2*3,lat0+size+0.05,'km',horizontalalignment = 'center')# style 1# print(help(ax.vlines))# ax.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=2, label='%d km' % (length))# ax.vlines(x = lon0, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=2)# ax.vlines(x = lon0+length/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=2)# # ax.text(lon0+length/2/111,lat0+size,'500 km',horizontalalignment = 'center')# ax.text(lon0+length/2/111,lat0+size,'%d' % (length/2),horizontalalignment = 'center')# ax.text(lon0,lat0+size,'0',horizontalalignment = 'center')# ax.text(lon0+length/111/2*3,lat0+size,'km',horizontalalignment = 'center')# style 2# plt.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=1, label='%d km' % (length))# plt.vlines(x = lon0, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)# plt.vlines(x = lon0+length/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)# plt.text(lon0+length/111,lat0+size,'%d km' % (length),horizontalalignment = 'center')# plt.text(lon0,lat0+size,'0',horizontalalignment = 'center')

添加比例尺

这个添加指北针的代码，原始代码是谁分享的无从考证！！！！

def add_north(ax, labelsize=18, loc_x=0.88, loc_y=0.85, width=0.06, height=0.09, pad=0.14):"""画一个比例尺带'N'文字注释主要参数如下:param ax: 要画的坐标区域 Axes实例 plt.gca()获取即可:param labelsize: 显示'N'文字的大小:param loc_x: 以文字下部为中心的占整个ax横向比例:param loc_y: 以文字下部为中心的占整个ax纵向比例:param width: 指南针占ax比例宽度:param height: 指南针占ax比例高度:param pad: 文字符号占ax比例间隙:return: None"""minx, maxx = ax.get_xlim()miny, maxy = ax.get_ylim()ylen = maxy - minyxlen = maxx - minxleft = [minx + xlen*(loc_x - width*.5), miny + ylen*(loc_y - pad)]right = [minx + xlen*(loc_x + width*.5), miny + ylen*(loc_y - pad)]top = [minx + xlen*loc_x, miny + ylen*(loc_y - pad + height)]center = [minx + xlen*loc_x, left[1] + (top[1] - left[1])*.4]triangle = mpatches.Polygon([left, top, right, center], color='k')ax.text(s='N',x=minx + xlen*loc_x,y=miny + ylen*(loc_y - pad + height),fontsize=labelsize,horizontalalignment='center',verticalalignment='bottom')ax.add_patch(triangle)

图像裁剪代码

这一段代码在咱们之前分享的博文中有涉及到！！具体可以看看这篇（点我）。

def shp2clip(originfig, ax, shpfile):'''originfig: colorbarax: 坐标轴shpfile: shp文件'''sf = shapefile.Reader(shpfile)vertices = []codes = []for shape_rec in sf.shapeRecords():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((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 originfig.collections:contour.set_clip_path(clip)return contour

栅格数据读取

数据用的还是咱们之前文章生成的fake江苏省气温数据，文末我给一个百度网盘链接吧，我把数据分享给大家，让大家好做测试！

库主要是使用的GDAL这个库，这个库如果大家有anaconda的话直接使用anaconda进行安装即可，如果没有的话，可以从这个网站上下载一下

values = gdal.Open('D:\CSDN\克里金插值/测试数据.tif')
x_ = values.RasterXSize  # 宽，读取一下x坐标有几个格点
y_ = values.RasterYSize  # 高，读取一下y坐标有几个格点
adfGeoTransform = values.GetGeoTransform() # 获取仿射矩阵
values = values.ReadAsArray() # 读取数据
# values_mask=np.ma.masked_where(values==0,values) #对0值进行mask
x = []
# 接下来这两个循环总体意思就是生成X,Y坐标（一维的！！）
for i in range(x_): x.append(adfGeoTransform[0] + i * adfGeoTransform[1]) #横坐标
y = []
for i in range(y_):y.append(adfGeoTransform[3] + i * adfGeoTransform[5]) #纵坐标
# print(adfGeoTransform)

画图函数

crs = ccrs.PlateCarree() # 设置投影
fig = plt.figure(figsize = (10, 15), dpi = 300) #创建一个绘图对象
ax1 = plt.subplot(1, 1, 1, projection = crs) #创建一个子图
geom = sr.Reader(r"D:\CSDN\克里金插值\江苏shp/江苏.shp").geometries() #读取shp文件
ax1.add_geometries(geom, crs,facecolor='none', edgecolor='black',linewidth=0.5) #绘制图形
ax1.add_feature(cfeature.OCEAN.with_scale('50m')) # 添加海洋
ax1.add_feature(cfeature.LAND.with_scale('50m')) # 添加陆地
ax1.add_feature(cfeature.RIVERS.with_scale('50m')) # 添加河流
ax1.add_feature(cfeature.LAKES.with_scale('50m')) # 添加湖泊
ax1.set_extent([116, 123, 30, 36]) # 设置显示范围
c = ax1.contourf(x, y, values, cmap='coolwarm',levels=np.arange(23, 28, 0.5),projection=crs) # 绘制等值线
gl = ax1.gridlines(draw_labels=True, linewidth=0.5, color='k', alpha=0.5, linestyle='--') # 设置网格线
# 如果不喜欢网格线，可以将上面的 linewidth=0.5 换成 linewidth=0
gl.xlabels_top = False
gl.ylabels_right = False
add_north(ax1)
add_scalebar(ax1,116.2,30.5,200,size=0.2) # 添加比例尺
shp2clip(c, ax1, r'D:\CSDN\克里金插值\江苏shp/江苏.shp') # 添加插值区域
plt.colorbar(c) # 添加颜色标尺
plt.show() # 显示图像

完整代码展示

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as cor
import cartopy.io.shapereader as sr
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import cartopy.crs as ccrs
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapefiledef add_north(ax, labelsize=18, loc_x=0.88, loc_y=0.85, width=0.06, height=0.09, pad=0.14):"""画一个比例尺带'N'文字注释主要参数如下:param ax: 要画的坐标区域 Axes实例 plt.gca()获取即可:param labelsize: 显示'N'文字的大小:param loc_x: 以文字下部为中心的占整个ax横向比例:param loc_y: 以文字下部为中心的占整个ax纵向比例:param width: 指南针占ax比例宽度:param height: 指南针占ax比例高度:param pad: 文字符号占ax比例间隙:return: None"""minx, maxx = ax.get_xlim()miny, maxy = ax.get_ylim()ylen = maxy - minyxlen = maxx - minxleft = [minx + xlen*(loc_x - width*.5), miny + ylen*(loc_y - pad)]right = [minx + xlen*(loc_x + width*.5), miny + ylen*(loc_y - pad)]top = [minx + xlen*loc_x, miny + ylen*(loc_y - pad + height)]center = [minx + xlen*loc_x, left[1] + (top[1] - left[1])*.4]triangle = mpatches.Polygon([left, top, right, center], color='k')ax.text(s='N',x=minx + xlen*loc_x,y=miny + ylen*(loc_y - pad + height),fontsize=labelsize,horizontalalignment='center',verticalalignment='bottom')ax.add_patch(triangle)#-----------函数：添加比例尺--------------
def add_scalebar(ax,lon0,lat0,length,size=0.45):'''ax: 坐标轴lon0: 经度lat0: 纬度length: 长度size: 控制粗细和距离的'''# style 3ax.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=1, label='%d km' % (length))ax.vlines(x = lon0, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)ax.vlines(x = lon0+length/2/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)ax.vlines(x = lon0+length/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)ax.text(lon0+length/111,lat0+size+0.05,'%d' % (length),horizontalalignment = 'center')ax.text(lon0+length/2/111,lat0+size+0.05,'%d' % (length/2),horizontalalignment = 'center')ax.text(lon0,lat0+size+0.05,'0',horizontalalignment = 'center')ax.text(lon0+length/111/2*3,lat0+size+0.05,'km',horizontalalignment = 'center')# style 1# print(help(ax.vlines))# ax.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=2, label='%d km' % (length))# ax.vlines(x = lon0, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=2)# ax.vlines(x = lon0+length/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=2)# # ax.text(lon0+length/2/111,lat0+size,'500 km',horizontalalignment = 'center')# ax.text(lon0+length/2/111,lat0+size,'%d' % (length/2),horizontalalignment = 'center')# ax.text(lon0,lat0+size,'0',horizontalalignment = 'center')# ax.text(lon0+length/111/2*3,lat0+size,'km',horizontalalignment = 'center')# style 2# plt.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=1, label='%d km' % (length))# plt.vlines(x = lon0, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)# plt.vlines(x = lon0+length/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1)# plt.text(lon0+length/111,lat0+size,'%d km' % (length),horizontalalignment = 'center')# plt.text(lon0,lat0+size,'0',horizontalalignment = 'center')def shp2clip(originfig, ax, shpfile):'''originfig: colorbarax: 坐标轴shpfile: shp文件'''sf = shapefile.Reader(shpfile)vertices = []codes = []for shape_rec in sf.shapeRecords():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((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 originfig.collections:contour.set_clip_path(clip)return contourvalues = gdal.Open('D:\CSDN\克里金插值/测试数据.tif')
x_ = values.RasterXSize  # 宽
y_ = values.RasterYSize  # 高
adfGeoTransform = values.GetGeoTransform() # 获取仿射矩阵values = values.ReadAsArray() # 读取数据
# values_mask=np.ma.masked_where(values==0,values) #对0值进行mask
x = []
for i in range(x_): x.append(adfGeoTransform[0] + i * adfGeoTransform[1]) #横坐标
y = []
for i in range(y_):y.append(adfGeoTransform[3] + i * adfGeoTransform[5]) #纵坐标
print(adfGeoTransform)crs = ccrs.PlateCarree()
fig = plt.figure(figsize = (10, 15), dpi = 300) #创建一个绘图对象
ax1 = plt.subplot(1, 1, 1, projection = crs) #创建一个子图
geom = sr.Reader(r"D:\CSDN\克里金插值\江苏shp/江苏.shp").geometries() #读取shp文件
ax1.add_geometries(geom, crs,facecolor='none', edgecolor='black',linewidth=0.5) #绘制图形
ax1.add_feature(cfeature.OCEAN.with_scale('50m')) # 添加海洋
ax1.add_feature(cfeature.LAND.with_scale('50m')) # 添加陆地
ax1.add_feature(cfeature.RIVERS.with_scale('50m')) # 添加河流
ax1.add_feature(cfeature.LAKES.with_scale('50m')) # 添加湖泊
ax1.set_extent([116, 123, 30, 36]) # 设置显示范围
c = ax1.contourf(x, y, values, cmap='coolwarm',levels=np.arange(23, 28, 0.5),projection=crs) # 绘制等值线
gl = ax1.gridlines(draw_labels=True, linewidth=0.5, color='k', alpha=0.5, linestyle='--') # 设置网格线
# 如果不喜欢网格线，可以将上面的 linewidth=0.5 换成 linewidth=0
gl.xlabels_top = False
gl.ylabels_right = False
add_north(ax1)
add_scalebar(ax1,116.2,30.5,200,size=0.2) # 添加比例尺
shp2clip(c, ax1, r'D:\CSDN\克里金插值\江苏shp/江苏.shp') # 添加插值区域
plt.colorbar(c) # 添加颜色标尺
plt.show() # 显示图像

总体出出来的图就是我下面这个样子的，如果大家不是很喜欢地图上的标记可以直接把上面添加湖泊、河流的那些代码注释掉就行了！

数据：
链接：https://pan.baidu.com/s/1uROH6QID2VswBl9Tu5nQ1w?pwd=GZWA
提取码：GZWA

博主说两句

上面出的图可能确实不完美，但是胜在是调用的开源库完成的，也不会涉及到任何版权问题。
如果大家还有什么建议的话，直接留言啦！
此外，最近有一个国产的气象数据处理库不错，官网在这里，这个库出图好看很多，但是我这边调用这个库会出来奇怪的错误！尚在研究之中。

用Python实现地理信息出图(含比例尺、指北针、图例)相关推荐

[ARCGIS]使用ARCPY+PYTHON+FME批量出图
应用场景: 最近在做某区退耕还林验收项目,要求按每个村出比例尺1:10000的JPG格式的图片,由于村比较多(接近400个,很多村在1:10000比例尺下无法显示完整,需要出多张图片),所以考虑批量完 ...
ArcMap出图小技巧：图例，比例尺，指北针，标题（附练习数据）
接上篇ArcMap出图小技巧分享讲解了主图的设计,接下来讲述一下出图时的其他要素(图例,比例尺,指北针,标题) 这时候我们要出图了,切换到布局视图. 系统默认纸张大小都是A4纸,如下图所示: 考虑到布 ...
【收藏】如何优雅的在 Python进行论文出图
python作为可视化利器,越来越受到研究人员的喜爱,并且常被作为论文出图的工具. 论文出图讲究固定的格式,所以配置一个固定的格式非常重要. 下面这篇博文列出了一些比较常见的一些配置和设置,供参考~ ...
python三维矩阵出图_python读取图片的方式,以及将图片以三维数组的形式输出方法...
python 三维npy数组如何画成三维图片画成三维图片? 你要是想要看空间分布的话画散点图就可以啊,用matplotlib 网页链接要是那种各种弯曲的面,也是matplotlib 网页链接再就 ...
使用Python辅助ArcGIS出图－使用ArcToolbox
目录 1.脚本编写 2.在ArcCatalog里面新建一个Toolbox 3.右键单击新建的Toolbox, add->script 4.选择脚本文件 5.设置输入和输出参数(选择待处理的文件夹 ...
arcgis 出图背景_ArcGIS空间制图分析视频教程（二狮兄出品）含ArcMap
这套教程是二狮兄出的一套ArcGIS地理空间制图数据分析视频教程,含ArcMap/ArcCatalog部分.教程分为上中下三部,已全部录制完毕,全部课程120节. 适用人群 ArcGIS目前的应用范围 ...
arcgis批量出图python代码_python使用arcpy.mapping模块批量出图
出图是项目里常见的任务,有的项目甚至会要上百张图片,所以批量出土工具很有必要.arcpy.mapping就是ArcGIS里的出图模块,能快速完成一个出图工具. arcpy.mapping模块里常用的类 ...
python画不出来图是什么原因-完美解决ARIMA模型中plot_acf画不出图的问题
问题描述:在画时间序列ACF时,调用 from statsmodels.graphics.tsaplots import plot_acf, plot_pacf plot_acf(data, lags ...
arcgis python脚本实现从界面选择输入输出_ArcGIS Python脚本实现数据驱动页面的批量出图...
这里讲一下如何通过Python脚本实现数据驱动页面的批量出图. 1 前提条件首先要知道的是数据驱动页面的导出必须启用驱动页面,同时地图文档必须处于布局视图中. 2 具体实现准备好批量出图的Pyth ...

用Python实现地理信息出图(含比例尺、指北针、图例)

文章目录

前言

库函数准备

分段讲解

添加比例尺

添加比例尺

图像裁剪代码

栅格数据读取

画图函数

完整代码展示

博主说两句

用Python实现地理信息出图(含比例尺、指北针、图例)相关推荐

最新文章

热门文章