jjzjj

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

卷心没有菜 2023-06-30 原文 
哈喽、哈喽大家好！！
最近浏览了不少代码，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

博主说两句

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

出图地理spanclasstokenpython开发语言人工智能

有关用Python实现地理信息出图(含比例尺、指北针、图例)的更多相关文章

  1. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声？ - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题？更新问题，以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点？我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等)，但我确实想创建一个输出文件。

  2. ruby-on-rails - Rails 常用字符串(用于通知和错误信息等) - 2

    大约一年前，我决定确保每个包含非唯一文本的Flash通知都将从模块中的方法中获取文本。我这样做的最初原因是为了避免一遍又一遍地输入相同的字符串。如果我想更改措辞，我可以在一个地方轻松完成，而且一遍又一遍地重复同一件事而出现拼写错误的可能性也会降低。我最终得到的是这样的:moduleMessagesdefformat_error_messages(errors)errors.map{|attribute,message|"Error:#{attribute.to_s.titleize}#{message}."}enddeferror_message_could_not_find(obje

  3. ruby - 解析 RDFa、微数据等的最佳方式是什么，使用统一的模式/词汇(例如 schema.org)存储和显示信息 - 2

    我主要使用Ruby来执行此操作，但到目前为止我的攻击计划如下:使用gemsrdf、rdf-rdfa和rdf-microdata或mida来解析给定任何URI的数据。我认为最好映射到像schema.org这样的统一模式，例如使用这个yaml文件，它试图描述数据词汇表和opengraph到schema.org之间的转换:#SchemaXtoschema.orgconversion#data-vocabularyDV:name:namestreet-address:streetAddressregion:addressRegionlocality:addressLocalityphoto:i

  4. ruby - 如何根据特征实现 FactoryGirl 的条件行为 - 2

    我有一个用户工厂。我希望默认情况下确认用户。但是鉴于unconfirmed特征，我不希望它们被确认。虽然我有一个基于实现细节而不是抽象的工作实现，但我想知道如何正确地做到这一点。factory:userdoafter(:create)do|user,evaluator|#unwantedimplementationdetailshereunlessFactoryGirl.factories[:user].defined_traits.map(&:name).include?(:unconfirmed)user.confirm!endendtrait:unconfirmeddoenden

  5. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意，但是谷歌搜索||=并不是很有帮助；)Python中是否有与Ruby和Perl中的||=语句等效的语句？例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外，类似这样的东西的通用术语是什么？条件分配是我的第一个猜测，但Wikipediapage跟我想的不太一样。

  6. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi？ - 2

    什么是ruby​​的rack或python的Java的wsgi？还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  7. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目：明明的随机数题目输入描述输出描述：示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od，od薪资待遇，od机试题清单华为OD机试真题大全，用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

  8. 【鸿蒙应用开发系列】- 获取系统设备信息以及版本API兼容调用方式 - 2

    在应用开发中，有时候我们需要获取系统的设备信息，用于数据上报和行为分析。那在鸿蒙系统中，我们应该怎么去获取设备的系统信息呢，比如说获取手机的系统版本号、手机的制造商、手机型号等数据。1、获取方式这里分为两种情况，一种是设备信息的获取，一种是系统信息的获取。1.1、获取设备信息获取设备信息，鸿蒙的SDK包为我们提供了DeviceInfo类，通过该类的一些静态方法，可以获取设备信息，DeviceInfo类的包路径为：ohos.system.DeviceInfo.具体的方法如下：ModifierandTypeMethodDescriptionstatic StringgetAbiList​()Obt

  9. python - 如何读取 MIDI 文件、更改其乐器并将其写回？ - 2

    我想解析一个已经存在的.mid文件，改变它的乐器，例如从“acousticgrandpiano”到“violin”，然后将它保存回去或作为另一个.mid文件。根据我在文档中看到的内容，该乐器通过program_change或patch_change指令进行了更改，但我找不到任何在已经存在的MIDI文件中执行此操作的库.他们似乎都只支持从头开始创建的MIDI文件。 最佳答案 MIDIpackage会为您完成此操作，但具体方法取决于midi文件的原始内容。一个MIDI文件由一个或多个音轨组成，每个音轨是十六个channel中任何一个上的

  10. 基于C#实现简易绘图工具【100010177】 - 2

    C#实现简易绘图工具一.引言实验目的:通过制作窗体应用程序(C#画图软件),熟悉基本的窗体设计过程以及控件设计,事件处理等,熟悉使用C#的winform窗体进行绘图的基本步骤,对于面向对象编程有更加深刻的体会.Tutorial任务设计一个具有基本功能的画图软件**·包括简单的新建文件,保存,重新绘图等功能**·实现一些基本图形的绘制,包括铅笔和基本形状等,学习橡皮工具的创建**·设计一个合理舒适的UI界面**注明:你可能需要先了解一些关于winform窗体应用程序绘图的基本知识,以及关于GDI+类和结构的知识二.实验环境Windows系统下的visualstudio2017C#窗体应用程序三.

随机推荐