«

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

时间:2024-8-5 09:34     作者:韩俊     分类: Python


本篇内容介绍了“Python如何实现批量读取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/h30v09.tif"
mcd_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/MCD12Q1.A2018001.h30v09.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即可。接下来的操作与本文开头提及的那两篇博客就一致,这里不再赘述~

标签: python

热门推荐