基於 Python 的地理空間繪圖(附源碼)

来源:https://www.cnblogs.com/1234567FENG/archive/2022/06/17/16386018.html
-Advertisement-
Play Games

前言 大部分情況下,地理繪圖可使用 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()

 

在這裡插入圖片描述

最後

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

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


您的分享是我們最大的動力!

-Advertisement-
Play Games
更多相關文章
  • 今天把平臺屬性的管理基本完成了,後臺管理做到現在基本也開始熟悉,確實就是對ElementUI的一個熟練程度。 一.平臺屬性管理 1.動態展示數據 先把介面弄好,應該在第三級標題選擇後進行發請求 靜態頁面搭建 渲染數據 屬性值列表,用到一個新組件 tag,並且這裡有多個屬性值,所以要遍歷,既然要在裡面 ...
  • 這裡給大家分享我在OpenLayers 地圖開發工作中總結出的一下代碼和註意點,希望對大家有所幫助 效果如下: 核心代碼展示:附帶講解註釋 var map = new ol.Map({ // 初始化地圖 target: 'map',// 選擇地圖對象 layers: [ new ol.layer.T ...
  • springboot整合ueditor實現圖片上傳和文件上傳功能 寫在前面: 在閱讀本篇之前,請先按照我的這篇隨筆完成對ueditor的前期配置工作: springboot+layui 整合百度富文本編輯器ueditor入門使用教程(踩過的坑)https://www.cnblogs.com/rain ...
  • 1、前言 作為嵌入式軟體開發,可能經常會使用命令行或者顯示屏等設備實現人機交互的功能,功能中通常情況都包含 UI 菜單設計;很多開發人員都會有自己的菜單框架模塊,防止重覆造輪子,網上有很多這種菜單框架的代碼,但是大多耦合性太強,無法獨立出來適配不同的菜單設計。 本文介紹一個降低了耦合性,完全獨立的菜 ...
  • 目錄 一.簡介 二.效果演示 三.源碼下載 四.猜你喜歡 零基礎 OpenGL (ES) 學習路線推薦 : OpenGL (ES) 學習目錄 >> OpenGL ES 基礎 零基礎 OpenGL (ES) 學習路線推薦 : OpenGL (ES) 學習目錄 >> OpenGL ES 轉場 零基礎 O ...
  • 分類方式 按參數分: 有參構造(預設構造) & 無參構造 按類型分: 普通構造 & 拷貝構造 調用方式 括弧法 顯示法 隱式轉換法 PS:下方所有文本均以此代碼為基礎 1 class Person { 2 public: 3 //無參構造函數 4 Person() { 5 std::cout << ...
  • 引言:今天工作遇到了一個需要按行讀取txt文件數據的需求,查詢了一下getline()函數,發現這竟然是一個C++的標準庫函數,而且設計的很好,特地做一下記錄。getline本質是一個定界流輸入截取函數,預設是換行符‘/n’ 個人技術博客(文章整理+源碼): https://zobolblog.gi ...
  • 大佬的理解->《IO流和File》 1、File類 File類是IO包中唯一代表磁碟文件本身的對象,File類定義了一些與平臺無關的方法來操作文件。通過調用File類提供的各種方法,能夠完成創建、刪除文件、重命名文件、判斷文件的讀寫許可權許可權是否存在、設置和查詢文件的最近修改時間等操作。 ​ File ...
一周排行
    -Advertisement-
    Play Games
  • 背景 在瀏覽器中訪問本地靜態資源html網頁時,可能會遇到跨域問題如圖。 是因為瀏覽器預設啟用了同源策略,即只允許載入與當前網頁具有相同源(協議、功能變數名稱和埠)的內容。 WebView2預設情況下啟用了瀏覽器的同源策略,即只允許載入與主機相同源的內容。所以如果我們把靜態資源發佈到iis或者通過node ...
  • 最近看幾個老項目的SQL條件中使用了1=1,想想自己也曾經這樣寫過,略有感觸,特別拿出來說道說道。編寫SQL語句就像炒菜,每一種調料的使用都會影響菜品的最終味道,每一個SQL條件的加入也會影響查詢的執行效率。那麼 1=1 存在什麼樣的問題呢?為什麼又會使用呢? ...
  • 好久不見,我又回來了。 給大家分享一個我最近使用c#代碼操作ftp伺服器的代碼示例: 1 public abstract class FtpOperation 2 { 3 /// <summary> 4 /// FTP伺服器地址 5 /// </summary> 6 private string f ...
  • 一:背景 1. 講故事 過年喝了不少酒,腦子不靈光了,停了將近一個月沒寫博客,今天就當新年開工寫一篇吧。 去年年初有位朋友找到我,說他們的系統會偶發性崩潰,在網上也發了不少帖子求助,沒找到自己滿意的答案,讓我看看有沒有什麼線索,看樣子這是一個牛皮蘚的問題,既然對方有了dump,那就分析起來吧。 二: ...
  • 自己製作的一個基於Entity Framework Core 的資料庫操作攔截器,可以列印資料庫執行sql,方便開發調試,代碼如下: /// <summary> /// EF Core 的資料庫操作攔截器,用於在資料庫操作過程中進行日誌記錄和監視。 /// </summary> /// <remar ...
  • 本文分享自華為雲社區《Go併發範式 流水線和優雅退出 Pipeline 與 Cancellation》,作者:張儉。 介紹 Go 的併發原語可以輕鬆構建流數據管道,從而高效利用 I/O 和多個 CPU。 本文展示了此類pipelines的示例,強調了操作失敗時出現的細微之處,並介紹了乾凈地處理失敗的 ...
  • 在上篇文章中,我們介紹到在多線程環境下,如果編程不當,可能會出現程式運行結果混亂的問題。出現這個原因主要是,JMM 中主記憶體和線程工作記憶體的數據不一致,以及多個線程執行時無序,共同導致的結果。 ...
  • 1、下載安裝包首先、進入官網下載安裝包網址:https://www.python.org/downloads/windows/下載步驟:進入下載地址,根據自己的電腦系統選擇相應的python版本 選擇適配64位操作系統的版本(查看自己的電腦操作系統版本), 點擊下載安裝包 也可以下載我百度雲分享的安 ...
  • 簡介 git-commit-id-maven-plugin 是一個maven 插件,用來在打包的時候將git-commit 信息打進jar中。 這樣做的好處是可以將發佈的某版本和對應的代碼關聯起來,方便查閱和線上項目的維護。至於它的作用,用官方說法,這個功能對於大型分散式項目來說是無價的。 功能 你 ...
  • 序言 在數字時代,圖像生成技術正日益成為人工智慧領域的熱點。 本討論將重點聚焦於兩個備受矚目的模型:DALL-E和其他主流AI繪圖方法。 我們將探討它們的優勢、局限性以及未來的發展方向。通過比較分析,我們期望能夠更全面地瞭解這些技術,為未來的研究和應用提供啟示。 Q: 介紹一下 dall-e Ope ...