MCD19A2 Python 处理

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 投影

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