专栏首页疯狂学习GISPython中gdal读取多波段HDF栅格影像并绘制直方图

Python中gdal读取多波段HDF栅格影像并绘制直方图

  本文详细介绍基于Python语言gdal等模块实现多波段HDF栅格图像文件(即.hdf文件)的读取、处理与像元值可视化等操作。此外,基于gdal等模块读取.tif格式栅格图层文件的方法可以查看Python中gdal实现多幅栅格影像批量绘制直方图,读取单波段.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_datasetIndex)是从0开始计算的;而后面的[0]则表示元组中的第一个参数,也就是上面一幅图中显示的该波段对应的数据路径。

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

欢迎关注公众号:疯狂学习GIS

本文分享自微信公众号 - 疯狂学习GIS(gh_c5ddfeeae7b4),作者:疯狂学习GIS

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2021-08-07

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 使用GDAL读取Sentinel数据

    版权声明:本文为博主原创文章,转载请注明原文出处! https://blog.csdn...

    卡尔曼和玻尔兹曼谁曼
  • GDAL读取MAIAC的HDF文件(熟肉)

    MAIAC全称Multi-Angle Implementationof Atmospheric Correction algorithm (多角度大气校正算法)...

    一个有趣的灵魂W
  • 读取HDF或者NetCDF格式的栅格数据

    HDF(Hierarchical Data Format)由NCSA(National Center for Supercomputing Applicatio...

    卡尔曼和玻尔兹曼谁曼
  • 栅格数据投影转换

    GDAL提供了gdalwarp命令可以方便地让我们进行影像拼接,重投影,裁剪,格式转换等功能

    卡尔曼和玻尔兹曼谁曼
  • GDAL简介

    Geospatial Data Abstraction Library (GDAL)是使用C/C++语言编写的用于读写空间数据的一套跨平台开源库。现有的大部分G...

    卡尔曼和玻尔兹曼谁曼
  • 使用Rasterio读取栅格数据

    有没有觉得用GDAL的Python绑定书写的代码很不Pythonic,强迫症的你可能有些忍受不了。不过,没关系,MapBox旗下的开源库Rasterio帮我们解...

    卡尔曼和玻尔兹曼谁曼
  • python+gdal+遥感图像拼接(mosaic)的实例

    作为摄影测量与遥感的从业者,笔者最近开始深入研究gdal,为工作打基础!个人觉得gdal也是没有什么技术含量,调用别人的api。但是想想这也是算法应用的一个技能...

    砸漏
  • python与地理空间分析(一)

    在气象数据分析中,地理空间要素是一个必须考虑的关键特征项,也是重要的影响因素。例如气温会随着海拔的升高而降低,地形的坡向朝向也会影响风速的分布,此外,典型的地形...

    zhangqibot
  • 打开栅格数据的正确方式

    下面的例子打开一副GeoTIFF影像,输出了影像的一些信息,然后遍历了所有波段,输出波段的一些信息

    卡尔曼和玻尔兹曼谁曼
  • 栅格数据创建与保存

    使用Python进行栅格数据处理,很多时候,我们会将GDAL的Dataset对象转化为NumPy的ndarray对象,这样我们可以使用很多通用的Python库对...

    卡尔曼和玻尔兹曼谁曼
  • 使用Python以优雅的方式实现根据shp数据对栅格影像进行切割

    一、前言        前面一篇文章(使用Python实现子区域数据分类统计)讲述了通过geopandas库实现对子区域数据的分类统计,说白了也就是如何根据一个...

    魏守峰
  • Maiac文件的AOD信息提取-by python

    Modis数据向来风骚,其HDF里包含了很多subdataset,其中有一个maiac的数据尤为特别。Maiac文件里含有大概12个数据集,每个数据集里又有4个...

    一个有趣的灵魂W
  • 栅格数据裁剪

    在进行遥感影像处理的时候,我们经常需要进行裁剪的工作,来看看如何使用GDAL工具进行这项操作吧!

    卡尔曼和玻尔兹曼谁曼
  • 栅格数据裁剪

    在进行遥感影像处理的时候,我们经常需要进行裁剪的工作,来看看如何使用GDAL工具进行这项操作吧!

    卡尔曼和玻尔兹曼谁曼
  • 如何使用Python处理HDF格式数据

    HDF也是一种自描述格式文件,主要用于存储和分发科学数据。气象领域中卫星数据经常使用此格式,比如MODIS,OMI,LIS/OTD等卫星产品。对HDF格式细节感...

    zhangqibot
  • 利用Python裁切tiff图像且读取tiff,shp文件的实例

    将转换成png后的图加载到软件中(专业软件ENVI5.3)查看结果详细信息如下图所示,成功的转换成png格式了。

    砸漏
  • 如何使用Python处理HDF格式数据及可视化问题

    原文链接:https://blog.csdn.net/Fairy_Nan/article/details/105914203

    砸漏
  • 如何使用Python处理HDF格式数据及可视化

    原文链接:https://blog.csdn.net/Fairy_Nan/article/details/105914203

    用户7886150
  • 使用Rasterio创建栅格数据

    在open()函数当中,我们可以像GDAL中的Create()方法一样,设置数据类型,数据尺寸,投影定义,仿射变换参数等一系列信息

    卡尔曼和玻尔兹曼谁曼

扫码关注云+社区

领取腾讯云代金券