基于 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級河流和湖泊為主，

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 = https://www.cnblogs.com/1234567FENG/archive/2022/06/17/[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 = https://www.cnblogs.com/1234567FENG/archive/2022/06/17/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()

最后

還有沒有學會的小伙伴嘛，點名批評不認真喲！關于今天的這一篇文章喜歡的記得點贊，讓我看看是哪一位靚仔在支持我，不懂

的也記得評論留言，學習的事馬虎不得，下一章見啦~

標籤：其他

上一篇：一元線性回歸的Python實作

下一篇：檔案輸入輸出處理-File