草庐IT

Python批量读取HDF多波段栅格数据并绘制像元直方图

fkxxgis 2023-03-28 原文

  本文介绍基于Python语言gdal模块，实现多波段HDF栅格图像文件的读取、处理与像元值可视化（直方图绘制）等操作。

  另外，基于gdal等模块读取.tif格式栅格图层文件的方法可以查看Python批量绘制遥感影像数据的直方图，读取单波段.hdf格式栅格图层文件的方法可以查看Python GDAL读取栅格数据并基于质量评估波段QA对指定数据加以筛选掩膜。

  本文期望实现的需求为：现有一存放.tif格式的全球LAI产品栅格数据的路径，需将这一路径下的全部LAI产品栅格数据依据另一路径下存放的全球MODIS植被覆盖类型产品栅格数据进行像元分类，并绘制全球每一种植被类型对应的LAI数值直方图。在这里，由于有前述两篇博客作为铺垫，本文对代码的讲解就着重于多波段HDF栅格图像文件的读取部分；其它内容由于在本文开头提及的前期两篇博客中已经详细介绍，这里就不再赘述~

  首先将本文所需代码展示如下：

# -*- coding: utf-8 -*-
"""
Created on Tue Jul 20 11:05:31 2021

@author: fkxxgis
"""

import os
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt

lai_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/h20v09.tif"
mcd_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/MCD12Q1.A2018001.h20v09.006.2019199233851.hdf"
pic_save_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/"

for veg_type in range(9):

    mcd_raster=gdal.Open(mcd_file_path)
    mcd_sub_dataset=mcd_raster.GetSubDatasets()
    hdf_band_num=len(mcd_sub_dataset)
    # for sub_dataset in mcd_sub_dataset:
    #     print(sub_dataset[1])
    # print(mcd_sub_dataset[2][1])
    mcd_sub_type=gdal.Open(mcd_sub_dataset[2][0])
    mcd_raster_array=mcd_sub_type.ReadAsArray()
    
    lai_raster=gdal.Open(lai_file_path)
    lai_raster_array=lai_raster.ReadAsArray()
    non_veg_type_lai_array=np.where(mcd_raster_array==veg_type+1,lai_raster_array,np.nan)
    plt.hist(non_veg_type_lai_array)
    plt.savefig(pic_save_path+"DRT_"+str(veg_type+1)+".png", dpi=300)
    plt.clf()
    plt.cla()

  我们直接讲解多波段HDF栅格图像文件读取部分的代码：首先，多波段.hdf格式文件的读取在一开始与单波段.hdf格式文件或.tif格式文件的读取一致，即通过gdal.Open()函数实现；但随后，需要额外借助len()函数获取HDF文件对应的波段数量。

  因为我们读取的HDF文件是多波段，因此hdf_band_num肯定是大于1的，那么刚刚读取进来的mcd_sub_dataset其实就是一个列表（List）；其中，这个列表的元素个数就是对应的多波段HDF文件波段数，列表的每一个元素则都是一个元组（tuple）；同时，每一个元组都有两个元素，其每一个元素都是一个字符串；其中第一个元素为当前HDF文件的当前波段对应的文件路径与部分提示信息，第二个元素作为当前HDF文件的当前波段对应的文件像素行列数、名称与数据类型。

  这么说可能不太明白，我们用一个实例来讲解。例如，通过上述代码读取一景具有六个波段的MODIS LAI产品——MCD15A3H产品，其第一个波段为FPAR数据，第二个波段为LAI数据。那么读取其后，得到的mcd_sub_dataset长这个样子：

  可以看到，是一个具有6个元素的列表。

  点开列表，可以看到6个元素每一个都是一个具有2个元素的元组：

  再点开第一个元组，可以看到其具有2个字符串格式的元素：

  其第二个元素包含了该波段对应的数据行数与列数（即[2400×2400]）、数据名称（即Fpar）、数据空间分辨率（即500m）、数据产品简称（即MOD_Grid_MCD15A3H），以及数据格式（即8-bit unsigned integer）；而第一个字符串没有显示完毕，我们可以点击打开看看：

  可以看到第一个元素则包含了该波段对应的数据路径、文件全称，以及部分与第二个元素重复的几个数据信息参数。

  有了上面的分析就比较清楚了，接下来再一次利用gdal.Open()函数读取我们需要的波段，mcd_sub_dataset[2][0]表示第三个波段；其中，第三个波段却用[2]来表示，是因为波段数量（也就是mcd_sub_dataset的Index）是从0开始计算的；而后面的[0]则表示元组中的第一个参数，也就是上面一幅图中显示的该波段对应的数据路径。

  随后，再利用.ReadAsArray()函数将其读取为Array即可。接下来的操作与本文开头提及的那两篇博客就一致，这里不再赘述~

栅格Pythoncodestrong

有关Python批量读取HDF多波段栅格数据并绘制像元直方图的更多相关文章

  1. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声？ - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题？更新问题，以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点？我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等)，但我确实想创建一个输出文件。

  2. ruby - 如何将脚本文件的末尾读取为数据文件(Perl 或任何其他语言) - 2

    我正在寻找执行以下操作的正确语法(在Perl、Shell或Ruby中):#variabletoaccessthedatalinesappendedasafileEND_OF_SCRIPT_MARKERrawdatastartshereanditcontinues. 最佳答案 Perl用__DATA__做这个:#!/usr/bin/perlusestrict;usewarnings;while(){print;}__DATA__Texttoprintgoeshere 关于ruby-如何将脚

  3. ruby - 解析 RDFa、微数据等的最佳方式是什么，使用统一的模式/词汇(例如 schema.org)存储和显示信息 - 2

    我主要使用Ruby来执行此操作，但到目前为止我的攻击计划如下:使用gemsrdf、rdf-rdfa和rdf-microdata或mida来解析给定任何URI的数据。我认为最好映射到像schema.org这样的统一模式，例如使用这个yaml文件，它试图描述数据词汇表和opengraph到schema.org之间的转换:#SchemaXtoschema.orgconversion#data-vocabularyDV:name:namestreet-address:streetAddressregion:addressRegionlocality:addressLocalityphoto:i

  4. Ruby 写入和读取对象到文件 - 2

    好的，所以我的目标是轻松地将一些数据保存到磁盘以备后用。您如何简单地写入然后读取一个对象？所以如果我有一个简单的类classCattr_accessor:a,:bdefinitialize(a,b)@a,@b=a,bendend所以如果我从中非常快地制作一个objobj=C.new("foo","bar")#justgaveitsomerandomvalues然后我可以把它变成一个kindaidstring=obj.to_s#whichreturns""我终于可以将此字符串打印到文件或其他内容中。我的问题是，我该如何再次将这个id变回一个对象？我知道我可以自己挑选信息并制作一个接受该信

  5. ruby - Ruby 有 `Pair` 数据类型吗？ - 2

    有时我需要处理键/值数据。我不喜欢使用数组，因为它们在大小上没有限制(很容易不小心添加超过2个项目，而且您最终需要稍后验证大小)。此外，0和1的索引变成了魔数(MagicNumber)，并且在传达含义方面做得很差(“当我说0时，我的意思是head...”)。散列也不合适，因为可能会不小心添加额外的条目。我写了下面的类来解决这个问题:classPairattr_accessor:head,:taildefinitialize(h,t)@head,@tail=h,tendend它工作得很好并且解决了问题，但我很想知道:Ruby标准库是否已经带有这样一个类？ 最佳

  6. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意，但是谷歌搜索||=并不是很有帮助；)Python中是否有与Ruby和Perl中的||=语句等效的语句？例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外，类似这样的东西的通用术语是什么？条件分配是我的第一个猜测，但Wikipediapage跟我想的不太一样。

  7. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi？ - 2

    什么是ruby​​的rack或python的Java的wsgi？还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  8. ruby - 我如何添加二进制数据来遏制 POST - 2

    我正在尝试使用Curbgem执行以下POST以解析云curl-XPOST\-H"X-Parse-Application-Id:PARSE_APP_ID"\-H"X-Parse-REST-API-Key:PARSE_API_KEY"\-H"Content-Type:image/jpeg"\--data-binary'@myPicture.jpg'\https://api.parse.com/1/files/pic.jpg用这个:curl=Curl::Easy.new("https://api.parse.com/1/files/lion.jpg")curl.multipart_form_

  9. 世界前沿3D开发引擎HOOPS全面讲解——集3D数据读取、3D图形渲染、3D数据发布于一体的全新3D应用开发工具 - 2

    无论您是想搭建桌面端、WEB端或者移动端APP应用，HOOPSPlatform组件都可以为您提供弹性的3D集成架构，同时，由工业领域3D技术专家组成的HOOPS技术团队也能为您提供技术支持服务。如果您的客户期望有一种在多个平台（桌面/WEB/APP，而且某些客户端是“瘦”客户端）快速、方便地将数据接入到3D应用系统的解决方案，并且当访问数据时，在各个平台上的性能和用户体验保持一致，HOOPSPlatform将帮助您完成。利用HOOPSPlatform，您可以开发在任何环境下的3D基础应用架构。HOOPSPlatform可以帮您打造3D创新型产品，HOOPSSDK包含的技术有：快速且准确的CAD

  10. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目：明明的随机数题目输入描述输出描述：示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od，od薪资待遇，od机试题清单华为OD机试真题大全，用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

随机推荐