前言 大部分情況下,地理繪圖可使用 Arcgis 等工具實現。但正版的 Arcgis 並非所有人可以承受。本文基於 Python 的 cartopy 和 matplotlib 等庫,為地理空間繪圖的代碼實現提供參考。 所有所需庫如下: gma、cartopy、matplotlib、numpy 更多內 ...
前言
大部分情況下,地理繪圖可使用 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()
最後
還有沒有學會的小伙伴嘛,點名批評不認真喲!關於今天的這一篇文章喜歡的記得點贊,讓我看看是哪一位靚仔在支持我,不懂
的也記得評論留言,學習的事馬虎不得。下一章見啦~