1. AI技术聚合首页
  2. Python

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

社会演员多 • • Python • 阅读 214 
哈喽、哈喽大家好！！
最近浏览了不少代码，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 3
    ax.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 - miny
    xlen = maxx - minx
    left = [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: colorbar
    ax: 坐标轴
    shpfile: shp文件
    '''
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        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((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 shapefile

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 - miny
    xlen = maxx - minx
    left = [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 3
    ax.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: colorbar
    ax: 坐标轴
    shpfile: shp文件
    '''
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        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((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

values = 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

博主说两句

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

文章出处登录后可见！

已经登录？立即刷新

共计人评分，平均分

到目前为止还没有投票！成为第一位评论此文章。

赞 (0)
社会演员多的头像社会演员多普通用户
0
上一篇 2023年4月5日
下一篇 2023年4月5日

相关推荐