Python实现批量读取HDF多波段栅格数据并绘制像元直方图
作者:疯狂学习GIS 发布时间:2023-08-20 17:27:08
本文介绍基于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即可。接下来的操作与本文开头提及的那两篇博客就一致,这里不再赘述~
来源:https://www.cnblogs.com/fkxxgis/p/17167617.html
猜你喜欢
- 前言一个表和多个表进行关联,但具体随着业务的加深,表不断的增加,关联的数量不断的增加,怎么通过一开始通过表的设计后,不在后期在修改表,彻底的
- 很多网友在浏览网页时应该会发现很多网页有显示时间和日期的功能,这个不难,使用可视化网页制作软件Drea
- 1、说明xlwt模块是非追加写入.xls模块,所以要一次性写入for循环和列表,这样就没有追加和非追加的说法。并且将Excel表合并,将每一
- 概述从今天开始我们将开启一段自然语言处理 (NLP) 的旅程. 自然语言处理可以让来处理, 理解, 以及运用人类的语言, 实现机器语言和人类
- 代码如下:--相信大家肯定经常会把数据导入到数据库中,但是可能会有些记录行的所有列的数据是null,这为null的数据是我们不需要 --现在
- 加了三个验证漏洞以及四个getshell方法# /usr/bin/env python3# -*- coding: utf-8 -*-# @
- python3.4以上的版本中,是默认自带pip的。查看pip的方法下载安装好python后,进入命令行,输入pip -V,即可查看pyth
- 1. 前言邮件,作为最正式规范的沟通方式,在日常办公过程中经常被用到我们都知道 Python内置了对 SMTP 的支持,可以发送纯文本、富文
- 前两天看的时候,所用的歌曲地址加密方式已变更。将以前的发出来供大家赏玩。解密函数是从flash里面反编译出来的,加密函数是自己根据解密函数写
- 步骤一:index页面处理<!DOCTYPE html><html lang="en"><
- 配置好virtualenv 和virtualenvwrapper后,使用pycharm创建新项目。之后要面临的问题就来了,之前一直使用的是s
- 人脸识别正在成为软件开发中的一种趋势。它有助于识别人脸并使应用程序更加健壮。在本教程中,我们将使用python和face_recogniti
- 我们在使用 requests 这类网络请求第三方库时,可以看到它有一个参数叫做 timeout ,就是指在网络请求发出开始计算,如果超过 t
- 如果你有一批IP地址想要获得这些IP具体的信息,比如归属国家,城市等,最好的办法当时是调用现有的api接口来获取,我在之前就写过一篇文章,是
- 写在前面的话:此篇还是asp相关的,相信玩ASP的都有这个感觉,当数据有5万多条时-------just like音乐网,要调用最新的10条
- 目录引言response响应元组形式响应make_response函数返回json格式数据其他特殊响应redirect() 重定向abort
- 一、概念梳理链表是计算机科学里面应用应用最广泛的数据结构之一。它是最简单的数据结构之一,同时也是比较高阶的数据结构(例如棧、环形缓冲和队列)
- Python使用Pika库(安装:sudo pip install pika)可以操作RabbitMQ消息队列服务器(安装:sudo apt
- 最近准备做一个关于scrapy框架的实战,爬取腾讯社招信息并存储,这篇博客记录一下创建项目的步骤pycharm是无法创建一个scrapy项目
- 想必大家都知道MSSQL中SA权限是什么,可以说是至高无上。今天我就它的危害再谈点儿,我所讲的是配合NBSI上传功能得到WebShell。在