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

哈喽、哈喽大家好!!
最近浏览了不少代码,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

博主说两句

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

  • 22
    点赞
  • 131
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 4
    评论
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

卷心没有菜

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值