MCD19A2 Python 处理

本文介绍了如何使用Python的gdal库处理MCD19A2遥感数据,包括导入库、读取数据、质量控制和保存处理结果。代码示例展示了读取AOD和AODQA数据,进行质量控制,并将HDF文件转换为GeoTIFF的过程。
摘要由CSDN通过智能技术生成

MCD19A2 Python 处理

# 处理MCD19A2数据需要一定的数据处理和编程能力,需要一定的学习和实践。如果您对Python基础没有很好的掌握,建议您先学习一些Python基础知识,例如变量、条件语句、循环、函数、列表和字典等,然后再逐步学习如何使用Python处理遥感数据。以下是处理MCD19A2数据的大致步骤:

# 1.导入Python库和模块,例如numpy、pandas和gdal等。
# 2.打开MCD19A2数据文件,使用gdal库读取数据集。
# 3.读取数据集中的时间、经度、纬度和AOD等数据,并将其保存到numpy数组或pandas数据框中。
# 4.对AOD数据进行质量控制,去除质量较差的数据,例如AODQA值低于某一阈值的数据。
# 5.将处理后的数据保存为文本文件或其他格式的数据文件,例如CSV、NetCDF等。
# 以下是一个处理MCD19A2数据的Python代码示例:

from osgeo import gdal
import numpy as np
import pandas as pd
import time

#将QA字段十进制转化为二进制  在进行位掩码运算
def convert_two(ten):
    '''数组十进制转二进制'''#转为二进制掩膜矩阵  最好质量为1 最差质量为0 进行掩膜
    two='{:0>16}'.format(str(bin(ten)[2:]))
    Best_quality=two[-12:-8]
    if Best_quality=='0000':
        result=int(1) 
    else:
        result=int(0)
    return result

# 打开MCD19A2数据文件
hdf_path=r'MCD19A2.A2023048.h33v11.006.2023050090446.hdf'
dataset = gdal.Open(hdf_path)#"MCD19A2.A2018001.h27v05.006.2018011223007.hdf"
# 读取AOD数据
SubDatasets=dataset.GetSubDatasets()#数据子集
# print(SubDatasets)
aod_band = SubDatasets[1][0]#grid1km:Optical_Depth_055 气溶胶550字符串
# print(aod_band)

#其中有多个波段 (4,1200,1200) 四是代表波段号
#以第一个波段为例
aod_dataset=gdal.Open(aod_band)
aod_data = np.array(aod_dataset.GetRasterBand(1).ReadAsArray())#1代表第一波段 从1开始 
# print(aod_data)

aodqa_band=SubDatasets[5][0]
# print(aodqa_band)
aodqa_dataset=gdal.Open(aodqa_band)
aodqa_data = np.array(aodqa_dataset.GetRasterBand(1).ReadAsArray())#1代表第一波段 从1开始 
# print(aodqa_data)

#对aod_data做质量控制
function_vector = np.vectorize(convert_two)#对矩阵做函数运算
mask_array=function_vector(aodqa_data)#转为二进制掩膜矩阵0-1矩阵
aod_qa_array=aod_data*mask_array#AOD*掩膜矩阵0-1矩阵
# print(aodqa_data)


#读取时间
Orbit_time_stamp=dataset.GetMetadata()['Orbit_time_stamp']#卫星白天过境时间
# print(Orbit_time_stamp)

Orbit_time_stamp_list=Orbit_time_stamp.split()#以空格分隔
# print(Orbit_time_stamp)

#读取波段1的时间
data_time=Orbit_time_stamp_list[0]
# print(data_time)

#判断是 Aqua Terra
if data_time.endswith('T'):#结尾是T 结尾
    print('Terra')
else:
    print('Aqua')


#时间读取转换
# 将字符串时间转换成时间组后在将其转换成时间戳格式
str_time=data_time[:-1]
# print(str_time) 
#'20230480210A' 
strptime=time.strptime(str_time,'%Y%j%H%M')#字符串转化为时间戳
# print(strptime)
# year=strptime.tm_year
# month=strptime.tm_mon
# mday=strptime.tm_mday
# hour=strptime.tm_hour
# mins=strptime.tm_min

#时间戳转为日期
time_list= time.strftime("%Y-%m-%d %H:%M:%S", strptime)
#print(time_list)

print(time_list, aod_data.flatten(),aodqa_data.flatten())



由于HDF文件里并没有直接给出经纬度信息, 但有相关的投影信息可以转换为sin 正弦投影,代码如下:

def HDF_to_GeoTiff(hdf_file_path,geotiff_outpath):
    dataset = gdal.Open(hdf_file_path)#打开文件
    subdatasets = dataset.GetSubDatasets()#获得子数据集
    str_Optical_Depth_055 = subdatasets[1][0]  # 气溶胶字符串
    aod_dataset=gdal.Open(str_Optical_Depth_055)#打开文件
    #aod_dataset = gdal.Open(path)
    im_width = aod_dataset.RasterXSize #栅格矩阵的列数
    im_height = aod_dataset.RasterYSize #栅格矩阵的行数
    im_bands = aod_dataset.RasterCount#栅格波段数
    im_geotrans=aod_dataset.GetGeoTransform()#获得仿射矩阵
    im_pro=aod_dataset.GetProjection()##获得空间投影
    driver = gdal.GetDriverByName('GTiff')
    #driver.CreateCopy(geotiff_outpath, aod_dataset)#, callback=progress, options=["TILED=YES", "COMPRESS=LZW"]
    dataset1=driver.Create(geotiff_outpath, im_width, im_height, im_bands, gdal.GDT_Int16)
    dataset1.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    dataset1.SetProjection(im_pro)  # 写入投影
    for num in range(im_bands):
        im_array=aod_dataset.GetRasterBand(num+1).ReadAsArray()
        im_array=im_array.astype(float)
        im_array[im_array<=0]=np.nan
        im_array[im_array>=5000]=np.nan
        im_array=np.nan_to_num(im_array)
        # print(im_array)
        dataset1.GetRasterBand(num+1).WriteArray(im_array)
        dataset1.GetRasterBand(num+1).SetNoDataValue(0)
        dataset1.GetRasterBand(num+1).FlushCache()
        dataset1.GetRasterBand(num+1).ComputeStatistics(False)#计算每个波段的统计数据 
        dataset1.BuildOverviews('average', [2,4,8,16,32])#创建概视图
    del dataset
    del aod_dataset
    del dataset1
    return geotiff_outpath +' ok !'

hdf_file_path=r'MCD19A2.A2023048.h33v11.006.2023050090446.hdf'
geotiff_outpath= r'MCD19A2.A2023048.h33v11.006.2023050090446_geo.tif'
HDF_to_GeoTiff(hdf_file_path,geotiff_outpath) #转化为sin 投影


希望以上能给大家带来帮助!

### 回答1: MCD19A2和MCTK是指MODIS(Moderate Resolution Imaging Spectroradiometer,中文为中分辨率成像光谱辐射计)卫星数据处理中的两个模块。 MCD19A2是MODIS的一个数据集,包含了大气辐射数据。它提供了地球大气层散射辐射和吸收辐射的信息。MCD19A2提供的数据可用于研究大气污染、气候变化、火灾等不同领域。使用MCD19A2数据可以进行大气辐射的监测和分析,为环境研究和防灾减灾提供帮助。 MCTK是MODIS的一个批处理工具包,用于对MODIS数据进行处理和分析。它包含了一系列的工具和算法,可以对MODIS的大气、陆地、海洋等不同类型的数据进行处理和计算。MCTK工具包中的批处理功能可以自动化处理大量的MODIS数据,提高处理效率和精度。 对于MCD19A2数据的批处理,可以使用MCTK工具包中的相关模块来完成。首先,需要进行数据的预处理,包括数据的下载、解压缩和格式转换等。接下来,可以使用MCTK提供的算法和工具对MCD19A2数据进行大气辐射的提取和分析,得到所需的结果和产品。最后,可以对处理后的数据进行可视化和报告生成,以便更好地理解和使用这些数据。 总之,MCD19A2和MCTK是MODIS卫星数据处理中的重要模块,它们提供了对大气辐射数据进行处理和分析的功能。通过使用MCTK批处理工具包,可以高效地处理大量的MODIS数据,并得到相应的结果和产品,为相关研究和应用提供支持。 ### 回答2: MCD19A2 MCTK批处理是指对MCD19A2和MCTK数据进行批量处理的过程。MCD19A2是MODIS火灾产品的一个子产品,包含了全球范围内的火点信息,而MCTK是火灾产品的一个工具包,用于对火点信息进行分类和处理。 对MCD19A2和MCTK进行批处理的目的是为了快速高效地处理大量的火点数据。在批处理过程中,可以使用自动化的方式对MCD19A2数据进行提取和整理,并且利用MCTK工具包对火点进行分类,例如将火点分为不同的类型(如林火、农田火等),并且计算火点强度等指标。 对火点数据进行批处理的好处是可以减轻人工处理的工作量,提高处理效率。通过批处理,可以快速地对大量的火点数据进行分析和统计,发现火灾的时空分布规律,为火灾监测和管理提供科学依据。 在进行MCD19A2 MCTK批处理时,需要先准备好MCD19A2和MCTK的数据,确保数据的完整性和准确性。然后,可以使用编程软件(如Python、MATLAB等)编写批处理代码,通过循环和条件判断来处理每一个火点数据,并将结果保存到指定的文件中。 总而言之,MCD19A2 MCTK批处理是对火点数据进行高效处理和分析的过程,能够大大提高工作效率和数据处理的准确性,为火灾监测和管理提供有力支持。
评论 14
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

香橙橙V

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值