Python中gdal实现MODIS卫星遥感影像产品栅格数据读取处理与质量控制QC波段筛选掩膜

简介: Python中gdal实现MODIS卫星遥感影像产品栅格数据读取处理与质量控制QC波段筛选掩膜

  下载后,打开HDF文件可以看到,其具有很多波段,同时包括质量控制QC波段;且在FPARLAI波段中,像元数值方面还具有精度较低的像元值、填充值等无效数值。上述这些都需要我们在读取数据时加以识别、处理与筛选。

  由于MODIS系列遥感影像产品种类较多,不同产品之间的属性差异较大;因此建议大家每次使用一种MODIS产品时,都到官网查看其基本信息,有需要的话还可以在官网下载对应产品的用户手册。前面提到,本文所用产品为MCD15A3H,因此可以在其官网https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/MCD15A3H#overview)查阅其基本信息或下载用户手册查看更为详细的产品属性。

  例如,下图所示即为用户手册中关于这一产品一景影像中波段分布情况与每一个波段具体信息的介绍表格;其中包括了波段含义、数据类型、填充值范围、有效值范围与缩放系数等关键参数,这些对于后期我们用gdal读取.hdf格式栅格文件而言具有重要意义。

  接下来,质量控制QC波段同样是执行栅格读取操作前有必要了解的信息。下图所示即为用户手册中关于这一产品一景影像中质量控制QC波段具体信息介绍的表格,其中包含了当前一景影像中FPARLAI产品的每一个像元所对应的算法、传感器、云覆盖等信息。这里需要注意的是:在MCD15A3H产品中是有两个质量控制QC波段的,这个是第一个QC,而第二个QC主要包括水陆区域、冰雪区域、气溶胶等信息,本文中暂且不涉及第二个QC。

  其中,由上表可知,QC波段的信息一共是由07共8个比特位(即Bit No.)组成,其中,由若干个比特位又可以组成Bit-word,每一个Bit-word就代表某一种QC波段信息。结合上图,我们可以对照下图这样一个实例进行理解:

  结合以上基本信息,我们已经对MCD15A3H产品的基本信息有了一定了解。接下来就可以进行栅格数据的读取与处理、筛选了。

  在这里需要注意的是,之前的两篇博客Python中gdal栅格遥感影像读取计算与写入处理及质量评估QA波段图层数据筛选掩膜https://blog.csdn.net/zhebushibiaoshifu/article/details/118878435)以及Python中gdal读取多波段HDF栅格遥感影像数据图层文件并依据像素绘制直方图
https://blog.csdn.net/zhebushibiaoshifu/article/details/119088429)已经对本次所要用到的大部分需求与代码加以实现并进行了详细讲解,这里就不再赘述。本文代码所实现功能与上述第一篇博客中的需求一致,唯一不同的是将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)

欢迎关注公众号/CSDN/知乎/微博:疯狂学习GIS


目录
打赏
0
1
1
0
116
分享
相关文章
【03】仿站技术之python技术,看完学会再也不用去购买收费工具了-修改整体页面做好安卓下载发给客户-并且开始提交网站公安备案-作为APP下载落地页文娱产品一定要备案-包括安卓android下载(简单)-ios苹果plist下载(稍微麻烦一丢丢)-优雅草卓伊凡
【03】仿站技术之python技术,看完学会再也不用去购买收费工具了-修改整体页面做好安卓下载发给客户-并且开始提交网站公安备案-作为APP下载落地页文娱产品一定要备案-包括安卓android下载(简单)-ios苹果plist下载(稍微麻烦一丢丢)-优雅草卓伊凡
42 13
【03】仿站技术之python技术,看完学会再也不用去购买收费工具了-修改整体页面做好安卓下载发给客户-并且开始提交网站公安备案-作为APP下载落地页文娱产品一定要备案-包括安卓android下载(简单)-ios苹果plist下载(稍微麻烦一丢丢)-优雅草卓伊凡
产品测评 | 上手分布式Python计算服务MaxFrame产品最佳实践
MaxFrame是阿里云自研的分布式计算框架,专为大数据处理设计,提供高效便捷的Python开发体验。其主要功能包括Python编程接口、直接利用MaxCompute资源、与MaxCompute Notebook集成及镜像管理功能。本文基于MaxFrame最佳实践,详细介绍了在DataWorks中使用MaxFrame创建数据源、PyODPS节点和MaxFrame会话的过程,并展示了如何通过MaxFrame实现分布式Pandas处理和大语言模型数据处理。测评反馈指出,虽然MaxFrame具备强大的数据处理能力,但在文档细节和新手友好性方面仍有改进空间。
🚀 MaxFrame 产品深度体验评测:Python 分布式计算的未来
在数据驱动的时代,大数据分析和AI模型训练对数据预处理的效率要求极高。传统的Pandas工具在小数据集下表现出色,但面对大规模数据时力不从心。阿里云推出的Python分布式计算框架MaxFrame,以“Pandas风格”为核心设计理念,旨在降低分布式计算门槛,同时支持超大规模数据处理。MaxFrame不仅保留了Pandas的操作习惯,还通过底层优化实现了高效的分布式调度、内存管理和容错机制,并深度集成阿里云大数据生态。本文将通过实践评测,全面解析MaxFrame的能力与价值,展示其在大数据和AI场景中的卓越表现。
63 4
🚀 MaxFrame 产品深度体验评测:Python 分布式计算的未来
MaxFrame 产品评测:大数据与AI融合的Python分布式计算框架
MaxFrame是阿里云MaxCompute推出的自研Python分布式计算框架,支持大规模数据处理与AI应用。它提供类似Pandas的API,简化开发流程,并兼容多种机器学习库,加速模型训练前的数据准备。MaxFrame融合大数据和AI,提升效率、促进协作、增强创新能力。尽管初次配置稍显复杂,但其强大的功能集、性能优化及开放性使其成为现代企业与研究机构的理想选择。未来有望进一步简化使用门槛并加强社区建设。
85 7
使用Python实现深度学习模型:智能质量检测与控制
使用Python实现深度学习模型:智能质量检测与控制 【10月更文挑战第8天】
442 62
使用Python实现深度学习模型:智能质量检测与控制
使用Python和PDFPlumber进行简历筛选:以SQL技能为例
本文介绍了一种使用Python和`pdfplumber`库自动筛选简历的方法,特别是针对包含“SQL”技能的简历。通过环境准备、代码解析等步骤,实现从指定文件夹中筛选出含有“SQL”关键词的简历,并将其移动到新的文件夹中,提高招聘效率。
75 8
使用Python和PDFPlumber进行简历筛选:以SQL技能为例
自动化测试的奥秘:如何用Selenium和Python提升软件质量
【9月更文挑战第35天】在软件开发的海洋中,自动化测试是那艘能引领我们穿越波涛的帆船。本文将揭开自动化测试的神秘面纱,以Selenium和Python为工具,展示如何构建一个简单而强大的自动化测试框架。我们将从基础出发,逐步深入到高级应用,让读者能够理解并实现自动化测试脚本,从而提升软件的质量与可靠性。
使用Python实现深度学习模型:智能药物研发与筛选
使用Python实现深度学习模型:智能药物研发与筛选
189 15
python django教学质量评价系统,实现学生、教师、管理员不同角色管理
本文介绍了一个基于Django框架开发的教学质量评价系统,该系统为学生、教师和管理员提供了不同角色的管理和评价功能,实现了教学质量的全方位评估和管理,旨在提高教育质量和促进教学改革。
126 5
python django教学质量评价系统,实现学生、教师、管理员不同角色管理
【2023年第十一届泰迪杯数据挖掘挑战赛】B题:产品订单的数据分析与需求预测 建模及python代码详解 问题一
本文详细介绍了第十一届泰迪杯数据挖掘挑战赛B题的解决方案,涵盖了对产品订单数据的深入分析、多种因素对需求量影响的探讨,并建立了数学模型进行未来需求量的预测,同时提供了Python代码实现和结果可视化的方法。
192 3
【2023年第十一届泰迪杯数据挖掘挑战赛】B题:产品订单的数据分析与需求预测 建模及python代码详解 问题一

热门文章

最新文章