基于 Python 的地理空间绘图详细教程

大部分情况下，地理绘图可使用 Arcgis 等工具实现。但正版的 Arcgis 并非所有人可以承受。本文基于 Python 的 cartopy 和 matplotlib 等库，为地理空间绘图的代码实现提供参考。

所有所需库如下：

gma、cartopy、matplotlib、numpy

Part1绘图目标

基于 Python 的地理空间绘图目标实现以下效果（包含比例尺、指北针、经纬网、图例等）：

Part2绘图思路

制图流程图

Part3数据处理

本例以 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)})

Part4绘制栅格

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)

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 = [ '林地', '灌木', '草地', '耕地', '建筑', '裸地/稀疏植被区', '雪和冰', '开阔水域',
'草本湿地', '红树林', '苔藓和地衣']

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配置制图范围

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]

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 = False
else:
gl.left_labels = False
if i < 2:
gl.bottom_labels = False
else:
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()

制图效果对比

文章转自：锐多宝的空间

点击关键词 ↓↓ 品往期精选

进阶手册！学习+考证

测绘师继续教育20学时

最后，欢迎来稿