前言
大部分情况下,地理绘图可使用 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级河流和湖泊为主。
python学习交流Q群:906715085### 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绘制数据
python学习交流Q群:906715085#### 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()
最后
还有没有学会的小伙伴嘛,点名批评不认真哟!关于今天的这一篇文章喜欢的记得点赞,让我看看是哪一位靓仔在支持我,不懂
的也记得评论留言,学习的事马虎不得。下一章见啦~
来源:https://www.cnblogs.com/1234567FENG/p/16386018.html
本站部分图文来源于网络,如有侵权请联系删除。