Python gdal讀取MODIS遙感影像並結合質量控制QC波段掩膜數據

来源:https://www.cnblogs.com/fkxxgis/archive/2023/03/10/17203063.html
-Advertisement-
Play Games

本文介紹基於Python中GDAL模塊，實現MODIS遙感影像數據的讀取、計算，並基於質量控制QC波段進行圖像掩膜的方法~ ...


  本文介紹基於Python中GDAL模塊，實現MODIS遙感影像數據的讀取、計算，並基於質量控制QC波段進行圖像掩膜的方法。

  前期的文章Python GDAL讀取柵格數據並基於質量評估波段QA對指定數據加以篩選掩膜詳細介紹了基於Python語言gdal等模塊實現遙感影像柵格數據的讀取，以及基於質量評估（QA）波段實現柵格像元篩選與掩膜的全部操作。而在本文，我們依據前述這一篇博客的代碼，結合大家更為熟悉的MODIS系列遙感影像產品，基於其質量評估波段進行具體的對照講解。也就是說，本文重點不在於代碼的講解（具體代碼在前述這一篇博客中已經很詳細地介紹了），而是將上述代碼在更為具體的一個實踐中加以應用，告訴大家該如何選擇波段、處理質量評估QA波段併進行篩選操作等。同時，這裡還有一點需要註意：在MODIS系列遙感影像中，質量評估波段更應該稱為質量控制波段，因為其官方手冊中將其寫作Quality Control，因此後文就寫作質量控制波段或QC波段。

  首先，需要下載好對應的MODIS數據，大家可以依據文章批量下載MODIS遙感影像：基於LAADS DAAC的方法中的方法進行下載。本文就以一景MODIS的LAI產品——MCD15A3H產品為例進行操作。

  下載後，打開HDF文件可以看到，其具有很多波段，同時包括質量控制QC波段；且在FPAR與LAI波段中，像元數值方面還具有精度較低的像元值、填充值等無效數值。上述這些都需要我們在讀取數據時加以識別、處理與篩選。

  由於MODIS系列遙感影像產品種類較多，不同產品之間的屬性差異較大；因此建議大家每次使用一種MODIS產品時，都到官網查看其基本信息，有需要的話還可以在官網下載對應產品的用戶手冊。前面提到，本文所用產品為MCD15A3H，因此可以在其官網查閱其基本信息或下載用戶手冊查看更為詳細的產品屬性。

  例如，下圖所示即為用戶手冊中關於這一產品一景影像中波段分佈情況與每一個波段具體信息的介紹表格；其中包括了波段含義、數據類型、填充值範圍、有效值範圍與縮放繫數等關鍵參數，這些對於後期我們用gdal讀取.hdf格式柵格文件而言具有重要意義。

  接下來，質量控制QC波段同樣是執行柵格讀取操作前有必要瞭解的信息。下圖所示即為用戶手冊中關於這一產品一景影像中質量控制QC波段具體信息介紹的表格，其中包含了當前一景影像中FPAR與LAI產品的每一個像元所對應的演算法、感測器、雲覆蓋等信息。這裡需要註意的是：在MCD15A3H產品中是有2個質量控制QC波段的，這個是第一個QC，而第二個QC主要包括水陸區域、冰雪區域、氣溶膠等信息，本文中暫且不涉及第二個QC。

  其中，由上表可知，QC波段的信息一共是由0至7共8個比特位（即Bit No.）組成，其中，由若幹個比特位又可以組成Bit-word，每一個Bit-word就代表某一種QC波段信息。結合上圖，我們可以對照下圖這樣一個實例進行理解：

  結合以上基本信息，我們已經對MCD15A3H產品的基本信息有了一定瞭解。接下來就可以進行柵格數據的讀取與處理、篩選了。

  在這裡需要註意的是，之前的兩篇文章Python GDAL讀取柵格數據並基於質量評估波段QA對指定數據加以篩選掩膜以及Python批量讀取HDF多波段柵格數據並繪製像元直方圖已經對本次所要用到的大部分需求與代碼加以實現併進行了詳細講解，這裡就不再贅述。本文代碼所實現功能與上述第一篇博客中的需求一致，唯一不同的是將GLASS產品更改為了MCD15A3H產品，且僅需對MCD15A3H產品的主演算法像元加以做差計算（也就是篩選出MCD15A3H產品中第一個QC波段對應二進位數的第一位為0的像元，其它像元就不用參與差值計算了）。

  具體代碼如下：

# -*- coding: utf-8 -*-
"""
Created on Sun Jul 25 14:57:45 2021

@author: fkxxgis
"""

import os
import copy
import numpy as np
from osgeo import gdal

rt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/RT_LAI/"
mcd15_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/mcd15A3H/"
out_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/"

rt_file_list=os.listdir(rt_file_path)
for rt_file in rt_file_list:
    rt_file_split=rt_file.split("_")
    rt_hv=rt_file_split[3][:-4]
    
    mcd15_file_list=os.listdir(mcd15_file_path)
    for mcd15_file in mcd15_file_list:
        if rt_hv in mcd15_file:
            rt_file_tif_path=rt_file_path+rt_file
            mcd15_file_tif_path=mcd15_file_path+mcd15_file
            
            drt_out_file_path=out_file_path+"drt/"
            if not os.path.exists(drt_out_file_path):
                os.makedirs(drt_out_file_path)
            drt_out_file_tif_path=drt_out_file_path+rt_hv+".tif"
            
            eco_out_file_path=out_file_path+"eco/"
            if not os.path.exists(eco_out_file_path):
                os.makedirs(eco_out_file_path)
            eco_out_file_tif_path=eco_out_file_path+rt_hv+".tif"
            
            wat_out_file_path=out_file_path+"wat/"
            if not os.path.exists(wat_out_file_path):
                os.makedirs(wat_out_file_path)
            wat_out_file_tif_path=wat_out_file_path+rt_hv+".tif"
            
            tim_out_file_path=out_file_path+"tim/"
            if not os.path.exists(tim_out_file_path):
                os.makedirs(tim_out_file_path)
            tim_out_file_tif_path=tim_out_file_path+rt_hv+".tif"
            
            rt_raster=gdal.Open(rt_file_tif_path)
            rt_raster_array=rt_raster.ReadAsArray()
            rt_lai_array=rt_raster_array[0]
            rt_qa_array=rt_raster_array[1]
            rt_lai_array_mask=np.where(rt_lai_array>30000,np.nan,rt_lai_array)
            rt_lai_array_fin=rt_lai_array_mask*0.001
            
            mcd15_raster=gdal.Open(mcd15_file_tif_path)
            mcd15_sub_dataset=mcd15_raster.GetSubDatasets()
            # for sub_dataset in mcd15_sub_dataset:
            #     print(sub_dataset[1])
            # print(mcd15_sub_dataset[1][1])
            # print(mcd15_sub_dataset[2][1])
            mcd15_sub_lai=gdal.Open(mcd15_sub_dataset[1][0])
            mcd15_sub_qc=gdal.Open(mcd15_sub_dataset[2][0])
            mcd15_lai_array=mcd15_sub_lai.ReadAsArray()
            mcd15_qc_array=mcd15_sub_qc.ReadAsArray()
            mcd15_lai_array_mask=np.where(mcd15_lai_array>248,np.nan,mcd15_lai_array)
            mcd15_lai_array_fin=mcd15_lai_array_mask*0.1
            
            rt_qa_array_bin=copy.copy(rt_qa_array)
            rt_qa_array_row,rt_qa_array_col=rt_qa_array.shape
            for i in range(rt_qa_array_row):
                for j in range(rt_qa_array_col):
                    rt_qa_array_bin[i][j]="{:012b}".format(rt_qa_array_bin[i][j])[-4:]
            
            mcd15_qc_array_bin=copy.copy(mcd15_qc_array)
            mcd15_qc_array_row,mcd15_qc_array_col=mcd15_qc_array.shape
            for i in range(mcd15_qc_array_row):
                for j in range(mcd15_qc_array_col):
                    mcd15_qc_array_bin[i][j]="{:08b}".format(mcd15_qc_array[i][j])[-1:]
            
            mcd15_lai_main_array=np.where(mcd15_qc_array_bin==1,np.nan,mcd15_lai_array_fin)
            
            lai_dif=rt_lai_array_fin-mcd15_lai_main_array
            lai_dif=lai_dif*1000
            
            drt_lai_dif_array=np.where((rt_qa_array_bin>=100) | (rt_qa_array_bin==11),
                                       np.nan,lai_dif)
            eco_lai_dif_array=np.where((rt_qa_array_bin<100) | (rt_qa_array_bin==111),
                                       np.nan,lai_dif)
            wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
                                       np.nan,lai_dif)
            tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
                                       np.nan,lai_dif)
            
            row=rt_raster.RasterYSize
            col=rt_raster.RasterXSize
            geotransform=rt_raster.GetGeoTransform()
            projection=rt_raster.GetProjection()
            
            # 輸出為int格式後，所得結果中0就是NoData
            driver=gdal.GetDriverByName("Gtiff")
            out_drt_lai=driver.Create(drt_out_file_tif_path,row,col,1,gdal.GDT_Int16)
            out_drt_lai.SetGeoTransform(geotransform)
            out_drt_lai.SetProjection(projection)
            out_drt_lai.GetRasterBand(1).WriteArray(drt_lai_dif_array)
            out_drt_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Int16)
            out_eco_lai.SetGeoTransform(geotransform)
            out_eco_lai.SetProjection(projection)
            out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
            out_eco_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Int16)
            out_wat_lai.SetGeoTransform(geotransform)
            out_wat_lai.SetProjection(projection)
            out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
            out_wat_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Int16)
            out_tim_lai.SetGeoTransform(geotransform)
            out_tim_lai.SetProjection(projection)
            out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
            out_tim_lai=None
            
            print(rt_hv)

  至此，大功告成。


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

-Advertisement-
Play Games
更多相關文章
  • python可變長參數
    哈嘍大家好 今天給大家分享一個用Python開發一款飛翔的小鳥游戲。 飛翔的小鳥（游戲英文名：Flappy Bird） 一款由越南獨立開發者開發的手機游戲，是之前非常流行的一款手機游戲 小游戲目標：讓小鳥穿過管子，不要碰到任何物體，挑戰更遠距離 今天，就讓我們一起用python來複刻一下這款游戲吧！ ...
  • SpringBoot 文件上傳+攔截器
    SpringBoot 文件上傳+攔截器 文件上傳原理 表單的enctype屬性規定在發送到伺服器之前應該如何對錶單數據進行編碼。 當表單的enctype="application/x-www-form-urlencoded"(預設)時，form表單中的數據格式為：key=value&key=valu ...
  • 推薦幾款主流好用的遠程終端連接管理軟體
    我可以明確地回答.我們之所以選擇Postgres，是因為它在操作上比MySQL更可靠，而當時公司的創始人相信SQL資料庫的可移植性. 隨著年份的發展，我們發現了這一點，我們發現基本上，Postgres是Rough中的這款鑽石，它具有一系列功能和一個開發社區，這是我們見過的最不可思議的開源項目之一，並 ...
  • Lua語言線上代碼運行編譯工具推薦
    Lua語言線上運行編譯,是一款可線上編程編輯器,在編輯器上輸入Lua語言代碼,點擊運行,可線上編譯運行Lua語言,Lua語言代碼線上運行調試,Lua語言線上編譯,可快速線上測試您的Lua語言代碼,線上編譯Lua語言代碼發現是否存在錯誤,如果代碼測試通過，將會輸出編譯後的結果。 該線上工具由IT寶庫提 ...
  • easyexcel填坑-自定義異常處理
    3.0版本以前listener中拋出的所有異常都會被easyexcel捕獲，並且包裝成ExcelAnalysisException，寫法如下： private void parseXmlSource(InputStream inputStream, ContentHandler handler) { ...
  • 如何通過C++ 將數據寫入 Excel 工作表
    直觀的界面、出色的計算功能和圖表工具，使Excel成為最流行的個人電腦數據處理軟體。在獨立的數據包含的信息量太少，而過多的數據又難以理清頭緒時，製作成表格是數據管理的最有效手段之一。這樣不僅可以方便整理數據，還可以方便我們查找和應用數據。後期我們還可以對具有相似表格框架，相同性質的數據進行合併彙總... ...
  • 快 40 歲，剛被裁。。
    作者：張飛洪 來源：www.cnblogs.com/jackyfei/p/13862607.html 離職的心態 人們在辭退或者被辭退都會對原公司抱有意見，因為疫情，公司業務告急，工資發不出來，我也失去了工作。雖然情緒上難免會有波動，但是轉念一想，我應該用開心的心態來看待這次辭職，並希望能快速翻過這 ...
  • 【機器學習】1. 廣義線性模型
    1. 廣義線性模型 1.1 從簡單線性回歸開始 機器學習的任務是從已知的數據中提取知識, 從而完成對新數據的識別與預測, 即應用知識. 但是, 數據本身不會給予研究者想要的答案, 大多數時候人們很難通過觀察或者簡單的計算提取有用的結論. 為了從歷史數據中整合知識, 人們提出了"模型"這一概念, 認為 ...
一周排行
    -Advertisement-
    Play Games
  • C#TMS系統學習（BaseCity頁面）
    C#TMS系統代碼-基礎頁面BaseCity學習 本人純新手，剛進公司跟領導報道，我說我是java全棧，他問我會不會C#，我說大學學過，他說這個TMS系統就給你來管了。外包已經把代碼給我了，這幾天先把增刪改查的代碼背一下，說不定後面就要趕鴨子上架了 Service頁面 //using => impo ...
  • C# - 委托與事件
    委托與事件 委托 委托的定義 委托是C#中的一種類型，用於存儲對方法的引用。它允許將方法作為參數傳遞給其他方法，實現回調、事件處理和動態調用等功能。通俗來講，就是委托包含方法的記憶體地址，方法匹配與委托相同的簽名，因此通過使用正確的參數類型來調用方法。 委托的特性 引用方法：委托允許存儲對方法的引用， ...
  • 解析 ABP vNext 依賴註入實現【屬性註入】的原理
    前言 這幾天閑來沒事看看ABP vNext的文檔和源碼，關於關於依賴註入(屬性註入)這塊兒產生了興趣。 我們都知道。Volo.ABP 依賴註入容器使用了第三方組件Autofac實現的。有三種註入方式，構造函數註入和方法註入和屬性註入。 ABP的屬性註入原則參考如下： 這時候我就開始疑惑了，因為我知道 ...
  • C#TMS系統學習（ShippingNotice頁面）
    C#TMS系統代碼-業務頁面ShippingNotice學習 學一個業務頁面，ok，領導開完會就被裁掉了，很突然啊，他收拾東西的時候我還以為他要旅游提前請假了，還在尋思為什麼回家連自己買的幾箱飲料都要叫跑腿帶走，怕被偷嗎？還好我在他開會之前拿了兩瓶芬達 感覺感覺前面的BaseCity差不太多，這邊的 ...
  • C#動態查詢：巧用Expression組合多條件表達式
    概述：在C#中，通過`Expression`類、`AndAlso`和`OrElse`方法可組合兩個`Expression<Func<T, bool>>`，實現多條件動態查詢。通過創建表達式樹，可輕鬆構建複雜的查詢條件。 在C#中，可以使用AndAlso和OrElse方法組合兩個Expression< ...
  • 在Biwen.QuickApi中整合一個極簡的發佈訂閱(事件匯流排)
    閑來無聊在我的Biwen.QuickApi中實現一下極簡的事件匯流排,其實代碼還是蠻簡單的,對於初學者可能有些幫助 就貼出來,有什麼不足的地方也歡迎板磚交流~ 首先定義一個事件約定的空介面 public interface IEvent{} 然後定義事件訂閱者介面 public interface I ...
  • WPF 穩定的全屏化視窗方法
    1. 案例 成某三甲醫預約系統， 該項目在2024年初進行上線測試，在正常運行了兩天後，業務系統報錯：The connection pool has been exhausted, either raise MaxPoolSize (currently 800) or Timeout (curren ...
  • 在 WPF 中集成 ASP.NET Core 和 WebView2 用於集成 SPA 應用
    背景 我們有些工具在 Web 版中已經有了很好的實踐，而在 WPF 中重新開發也是一種費時費力的操作，那麼直接集成則是最省事省力的方法了。 思路解釋 為什麼要使用 WPF？莫問為什麼，老 C# 開發的堅持，另外因為 Windows 上已經裝了 Webview2/edge 整體打包比 electron ...
  • C#事件
    EDP是一套集組織架構，許可權框架【功能許可權，操作許可權，數據訪問許可權，WebApi許可權】，自動化日誌，動態Interface，WebApi管理等基礎功能於一體的，基於.net的企業應用開發框架。通過友好的編碼方式實現數據行、列許可權的管控。 ...
  • .Net8.0 Blazor Hybird 桌面端 (WPF/Winform) 發佈到 Win7+
    .Net8.0 Blazor Hybird 桌面端 (WPF/Winform) 實測可以完整運行在 win7sp1/win10/win11. 如果用其他工具打包,還可以運行在mac/linux下, 傳送門BlazorHybrid 發佈為無依賴包方式 安裝 WebView2Runtime 1.57 M ...
所有分類