基于 Python 的地理空间绘图指南
大部分情况下,地理绘图可使用 Arcgis 等工具实现。但正版的 Arcgis 并非所有人可以承受。本文基于 Python 的 cartopy 和 matplotlib 等库,为地理空间绘图的代码实现提供参考。
所有所需库如下:
gma、cartopy、matplotlib、numpy
更多内容可转到:地理与气象分析库----使用指南。
0 绘图目标
基于 Python 的地理空间绘图目标实现以下效果(包含比例尺、指北针、经纬网、图例等):
1 绘图思路
2 数据处理
本例以 ESA 2020 陆表覆盖河南省地物分类数据为例,通过gma.rasp.AddColorTable 更新色彩映射表,形成三个与原始文件不同的副本栅格(仅配色不同)。并对四个栅格进行绘制。这四个文件分别为:
“地表覆盖_河南_ESA_2020.tif” ----原始数据
“地表覆盖_河南_ESA_2020 - 副本.tif”
“地表覆盖_河南_ESA_2020 - 副本 (2).tif”
“地表覆盖_河南_ESA_2020 - 副本 (3).tif”
底图以我国省、地级市边界以及1-5级河流和湖泊为主。
import gma
# 1.根据定义更新——第一个副本
## 待更新的色彩映射表
ColorTable = {10:(0,112,255,255),20:(255,211,127,255),30:(76,230,0,255),40:(123,104,238,255),50:(230,230,0,255),60:(205,245,122,255),70:(156,200,121,255),80:(245,162,122,255),90:(190,210,255,255),95:(109,150,178,255),100:(223,198,142,255)}
## 将定义的色彩映射表更新到 副本
gma.rasp.AddColorTable("地表覆盖_河南_ESA_2020 - 副本.tif",ColorTable = ColorTable)
# 2.根据模板栅格更新——第二个副本
## 将 副本 的色彩映射表更新到 副本(2)
gma.rasp.AddColorTable("地表覆盖_河南_ESA_2020 - 副本 (2).tif","地表覆盖_河南_ESA_2020 - 副本.tif")
# 3.根据模板栅格和定义更新——第三个副本
## 将 副本 以及定义的色彩映射表更新到 副本 (3)
gma.rasp.AddColorTable("地表覆盖_河南_ESA_2020 - 副本 (3).tif","地表覆盖_河南_ESA_2020 - 副本.tif",ColorTable = {10:(100,100,100,255), 40:(200,200,200,255)})
3 绘制栅格
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.colors as cor
import numpy as np
import gma.extend.mapplottools as mpt
PAR = {'font.sans-serif': 'Times New Roman','axes.unicode_minus': False,}
plt.rcParams.update(PAR)
3.1 读取色彩映射表信息(若不包含,可自行定义色带)
InFiles = ["地表覆盖_河南_ESA_2020.tif", "地表覆盖_河南_ESA_2020 - 副本.tif", "地表覆盖_河南_ESA_2020 - 副本 (2).tif", "地表覆盖_河南_ESA_2020 - 副本 (3).tif"]
#### 读取四组数据色彩信息
CMap = []
Colors = []
for InFile in InFiles:DataSet = gma.Open(InFile)Color = DataSet.GetGDALDataset().GetRasterBand(1).GetColorTable()ColorTable = [Color.GetColorEntry(i) for i in range(Color.GetCount())]# 转换 色彩映射表 为 Matplotlib 可识别的格式CMapV = tuple(tuple(np.array(CT) / 255) for CT in ColorTable)# 生成色带CMap.append(cor.ListedColormap(CMapV))Colors.append([CMapV[i] for i in range(10, 110, 10)] + [CMapV[95]])
#### 为四组数据分配名称
Method = ['原始配色', '根据定义更新', '根据模板栅格更新', '根据模板栅格和定义更新']
#### 为颜色值定义含义
ColorName = ['林地', '灌木', '草地', '耕地', '建筑', '裸地/稀疏植被区', '雪和冰', '开阔水域', '草本湿地', '红树林', '苔藓和地衣']
3.2 定义数据分块——用于数据分块绘制(节约内存)
当数据过大时,直接绘制可能失败。若想精确绘制,可采用此方法(若涉及到投影,大数据耗时较久)。否则,可以缩放数据,减小分辨率(类似栅格金字塔构建规则)进行绘制。
BlockSize = 8000
Columns = DataSet.Columns
Rows = DataSet.Rows
Blocks = [(r, c) for r in range(0, Rows, BlockSize) for c in range(0, Columns, BlockSize)]
3.3 配置制图范围
GEOT = DataSet.GeoTransform
Columns = DataSet.Columns
Rows = DataSet.Rows
# 数据边界
ExtentData = [GEOT[0], GEOT[0] + GEOT[1] * Columns, GEOT[3] + GEOT[-1] * Rows, GEOT[3]]
# 绘图边界(以数据边界为基础确定)
EL, ER, EB, ET = 0.2, 0.1, 0.15, 0.05 # 左右、下上边界的扩展比例
ExtentPLT = [ExtentData[0] - (ExtentData[1] - ExtentData[0]) * EL, ExtentData[1] + (ExtentData[1] - ExtentData[0]) * ER, ExtentData[2] - (ExtentData[3] - ExtentData[2]) * EB, ExtentData[3] + (ExtentData[3] - ExtentData[2]) * ET]
3.4 绘制数据
WKTCRS = DataSet.Projection
DataCRS = mpt.GetCRS(WKTCRS)
fig = plt.figure(figsize = (10, 10), dpi = 600)# 定义一个标准中国区 ALBERS 投影
Alberts_China = ccrs.AlbersEqualArea(central_longitude = 105, standard_parallels = (25.0, 47.0)) for i in range(4):ax = plt.subplot(2, 2, i + 1, projection = Alberts_China) # 0.控制数据显示范围ax.set_extent(ExtentPLT, crs = DataCRS)# 1.绘制底图图层(应用自有高精度数据做底图)## 1.1 添加行政边界mpt.AddGeometries(ax, r"Region\VTD_PG_PLCity_China.shp", EdgeColor = 'LightGrey', LineWidth = 0.1)mpt.AddGeometries(ax, r"Region\VTD_PG_Province_China.shp", EdgeColor = 'Gray', LineWidth = 0.2)## 1.2 添加河流湖泊mpt.AddGeometries(ax, r"river\1级河流.shp", EdgeColor = 'RoyalBlue', LineWidth = 0.4)mpt.AddGeometries(ax, r"river\2级河流.shp", EdgeColor = 'DodgerBlue', LineWidth = 0.3)mpt.AddGeometries(ax, r"river\3级河流.shp", EdgeColor = 'DeepSkyBlue', LineWidth = 0.2)mpt.AddGeometries(ax, r"river\4级河流.shp", EdgeColor = 'SkyBlue', LineWidth = 0.15)mpt.AddGeometries(ax, r"river\5级河流.shp", EdgeColor = 'LightSkyBlue', LineWidth = 0.05)mpt.AddGeometries(ax, r"river\主要湖泊.shp", EdgeColor = 'none', LineWidth = 0, FaceColor = '#BEE8FF')# 2.绘制数据图层## 分块绘制(节约内存)for Block in Blocks:# 数据都一样,读取一个文件的数据即可DrawData = DataSet.ToArray(*Block, BlockSize, BlockSize)ExtentBlock = [GEOT[0] + Block[1] * GEOT[1], GEOT[0] + (DrawData.shape[1] + Block[1]) * GEOT[1], GEOT[3] - (DrawData.shape[0] + Block[0]) * GEOT[1], GEOT[3] - Block[0] * GEOT[1]]im = ax.imshow(DrawData, transform = DataCRS, cmap = CMap[i], extent = ExtentBlock, zorder = 2,interpolation = 'none', vmin = 0, vmax = 255)# 3.为绘制区域增加经纬网gl = ax.gridlines(draw_labels = True, dms = False, x_inline = False, y_inline = False, linestyle = (0, (10, 10)), linewidth = 0.2,color = 'Gray',rotate_labels = False,xlabel_style = {'fontsize': 8},ylabel_style = {'fontsize': 8})## 3.1忽略相邻轴的经纬网标签if i % 2 == 0:gl.right_labels = Falseelse:gl.left_labels = Falseif i < 2:gl.bottom_labels = Falseelse:gl.top_labels = False ax.set_title(Method[i], fontsize = 10, y = 0.92, fontdict = {'family':'SimSun'})# n.其他优化设置## n.1 添加指北针mpt.AddCompass(ax, LOC = (0.2, 0.85), SCA = 0.04, FontSize = 10)## n.2 添加比例尺mpt.AddScaleBar(ax, LOC = (0.8, 0.08), SCA = 0.1, FontSize = 6, PROJType = 'PROJCS', UnitPad = 0.25, BarWidth = 0.6)## n.3 添加图例并修饰mpt.AddLegend(ax, Colors[i], LegendName = '分类', LengedInterval = 0.4, LabelList = ColorName, LegendSize = 8, TextInterval = 0.1, LOC = (0.05, 0.32), SCA = 0.03, AspectRatio = 1.5, Columns = 2, ColumnWide = 0.15, RowInterval = 0.015, FontSize = 6, EdgeColor = 'k', EdgeWidth = 0.1)
plt.subplots_adjust(wspace = 0.05, hspace = -0.05)
plt.show()
4 支持与赞助
基于 Python 的地理空间绘图指南相关推荐
- 【教程】基于 Python 的地理空间绘图指南
大部分情况下,地理绘图可使用 ArcGIS 等工具实现.但正版的 ArcGIS 并非所有人可以承受.本文基于 Python 的 cartopy 和 matplotlib 等库,为地理空间绘图的代码实现 ...
- python空间分析_读书笔记——《python地理空间分析指南》
本文为<Python地理空间分析指南(第2版)>的读书摘录,顺便挖个坑,进一步对python的几个包做学习整理. 本笔记的用途:了解python地理空间处理的技术框架和实现途径. 第三章 ...
- 《Python地理空间分析指南 第2版》学习笔记-5.1 距离测量
第5章 Python与地理信息系统 本章主要学习Python处理矢量数据,包含以下内容: 距离测量 坐标转换 矢量数据重投影 Shapefile 文件编辑 海量数据过滤 专题地图创建 非GIS数据类型 ...
- python地理空间分析指南pdf邓世超_Python地理空间分析指南(第2版)源代码.zip
[实例简介] Python地理空间分析指南(第2版)的随书源代码,需要的朋友可以下载一下~~ [实例截图] [核心代码] Python地理空间分析指南(第2版)源代码 └── Python地理空间分析 ...
- Python地理空间分析指南(第2版)学习笔记01
目录 前言 一.任务 二.实现与解析 1.引入库 2.构造数据模型 3.渲染地图元素 4.执行查询操作以及完成绘图 三.总结 前言 本书假定读者了解Pyhon.信息技术的基本知识,并且至少对地理空间分 ...
- 《Python地理空间分析指南(第2版)》——1.9 地理信息系统基本概念
本节书摘来自异步社区<Python地理空间分析指南(第2版)>一书中的第1章,第1.9节,作者: [美]Joel Lawhead(莱哈德) 更多章节内容可以访问云栖社区"异步社区 ...
- 一直在构建工作空间_国际资讯Python与地理空间分析
点击图片上方蓝色字体"慧天地"即可订阅 英文原文来源:www.gislounge.com 英文原文链接:https://www.gislounge.com/python-and-g ...
- gma 地理空间绘图:(1)绘制简单的世界地图-1.地图绘制与细节调整
了解 gma gma 是什么? gma 是一个基于 Python 的地理.气象数据快速处理和数据分析函数包(Geographic and Meteorological Analysis,gma).gm ...
- python空间数据处理_基于Python语言的空间数据处理
龙源期刊网 http://www.doczj.com/doc/7b0e0476172ded630a1cb662.html 基于Python语言的空间数据处理 作者:何丽娴甘淑陈应跃 来源:<价值 ...
最新文章
- 隔离太无聊?每天一个数据科学项目,数据集都准备好了!
- SAP独门神器之VC变式配置,硬核整理版重磅推出!
- 关于Verilog HDL的一些技巧、易错、易忘点(不定期更新)
- java soa例子_哪位大牛能举个实例讲下SOA与传统架构的区别?
- 异常检测——无监督、高斯分布模型,需要带标记的样本数据,基本假设:特征符合高斯分布...
- Java SSM4——Spring
- ArrayList详细
- 如何配置RadASM
- 使用Excel2016求解运筹学线性规划
- 手机隐藏Magisk的root痕迹,适用于含zygisk的Magisk
- Linux一句话精彩问答
- VsCode文件屏蔽
- /var/tmp/rpm-tmp. 安装失败时找不到tmp文件的应对方法
- 渠道商用假流量冒充真实用户
- ARM7——LPC2xxx小总结
- codeforces 596E Wilbur and Strings(DFS)
- 微信小程序canvas2d使用封装与案例使用
- 陕西省职业计算机考试试题,2010陕西省计算机等级考试试题 二级C试题最新考试试题库...
- The thirteen day
- vue使用jointJs,vue流程图、旅程图