基于 Python 的地理空间绘图指南
生活随笔
收集整理的這篇文章主要介紹了
基于 Python 的地理空间绘图指南
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
??大部分情況下,地理繪圖可使用 Arcgis 等工具實現。但正版的 Arcgis 并非所有人可以承受。本文基于 Python 的 cartopy 和 matplotlib 等庫,為地理空間繪圖的代碼實現提供參考。
??所有所需庫如下:
gma、cartopy、matplotlib、numpy
??更多內容可轉到:地理與氣象分析庫----使用指南。
0 繪圖目標
??基于 Python 的地理空間繪圖目標實現以下效果(包含比例尺、指北針、經緯網、圖例等):
1 繪圖思路
#mermaid-svg-2bTYQcEsAKsH3kRC {font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:16px;fill:#333;}#mermaid-svg-2bTYQcEsAKsH3kRC .error-icon{fill:#552222;}#mermaid-svg-2bTYQcEsAKsH3kRC .error-text{fill:#552222;stroke:#552222;}#mermaid-svg-2bTYQcEsAKsH3kRC .edge-thickness-normal{stroke-width:2px;}#mermaid-svg-2bTYQcEsAKsH3kRC .edge-thickness-thick{stroke-width:3.5px;}#mermaid-svg-2bTYQcEsAKsH3kRC .edge-pattern-solid{stroke-dasharray:0;}#mermaid-svg-2bTYQcEsAKsH3kRC .edge-pattern-dashed{stroke-dasharray:3;}#mermaid-svg-2bTYQcEsAKsH3kRC .edge-pattern-dotted{stroke-dasharray:2;}#mermaid-svg-2bTYQcEsAKsH3kRC .marker{fill:#333333;stroke:#333333;}#mermaid-svg-2bTYQcEsAKsH3kRC .marker.cross{stroke:#333333;}#mermaid-svg-2bTYQcEsAKsH3kRC svg{font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:16px;}#mermaid-svg-2bTYQcEsAKsH3kRC .label{font-family:"trebuchet ms",verdana,arial,sans-serif;color:#333;}#mermaid-svg-2bTYQcEsAKsH3kRC .cluster-label text{fill:#333;}#mermaid-svg-2bTYQcEsAKsH3kRC .cluster-label span{color:#333;}#mermaid-svg-2bTYQcEsAKsH3kRC .label text,#mermaid-svg-2bTYQcEsAKsH3kRC span{fill:#333;color:#333;}#mermaid-svg-2bTYQcEsAKsH3kRC .node rect,#mermaid-svg-2bTYQcEsAKsH3kRC .node circle,#mermaid-svg-2bTYQcEsAKsH3kRC .node ellipse,#mermaid-svg-2bTYQcEsAKsH3kRC .node polygon,#mermaid-svg-2bTYQcEsAKsH3kRC .node path{fill:#ECECFF;stroke:#9370DB;stroke-width:1px;}#mermaid-svg-2bTYQcEsAKsH3kRC .node .label{text-align:center;}#mermaid-svg-2bTYQcEsAKsH3kRC .node.clickable{cursor:pointer;}#mermaid-svg-2bTYQcEsAKsH3kRC .arrowheadPath{fill:#333333;}#mermaid-svg-2bTYQcEsAKsH3kRC .edgePath .path{stroke:#333333;stroke-width:2.0px;}#mermaid-svg-2bTYQcEsAKsH3kRC .flowchart-link{stroke:#333333;fill:none;}#mermaid-svg-2bTYQcEsAKsH3kRC .edgeLabel{background-color:#e8e8e8;text-align:center;}#mermaid-svg-2bTYQcEsAKsH3kRC .edgeLabel rect{opacity:0.5;background-color:#e8e8e8;fill:#e8e8e8;}#mermaid-svg-2bTYQcEsAKsH3kRC .cluster rect{fill:#ffffde;stroke:#aaaa33;stroke-width:1px;}#mermaid-svg-2bTYQcEsAKsH3kRC .cluster text{fill:#333;}#mermaid-svg-2bTYQcEsAKsH3kRC .cluster span{color:#333;}#mermaid-svg-2bTYQcEsAKsH3kRC div.mermaidTooltip{position:absolute;text-align:center;max-width:200px;padding:2px;font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:12px;background:hsl(80, 100%, 96.2745098039%);border:1px solid #aaaa33;border-radius:2px;pointer-events:none;z-index:100;}#mermaid-svg-2bTYQcEsAKsH3kRC :root{--mermaid-font-family:"trebuchet ms",verdana,arial,sans-serif;} 柵格 讀取柵格數據 底圖 讀取矢量底圖 確定繪圖區域 繪制數據 繪制底圖 添加指北針\比例尺\圖例 輸出繪圖結果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价值算法_第十课-Pytho