前言
上一篇文章有小伙伴问我，基于 AERONET 站点观测数据对比验证AOD数据产品,能不能利用Python代码来实现，以下提供以下我利用地基观测数据验证卫星数据的思路，本文以验证MYD04_L2 数据为例，供大家参考。
思路：遍历每一条地基观测数据，来搜寻对应时空窗口的匹配文件
时空匹配策略（见上一篇文章）：
1.地面观测点时间±30 min 之内的卫星观测值
2.地面观测站点 3*3像元平均值（MYD04_L2 分辨率约10km）
一、整理数据文件夹
为了方便对卫星过境时间的迅速匹配，需要先对数据文件夹进行整理，将其归并到属于自己的 …/年/日 文件目录下。
整理之前
移动文件整理代码
#!/usr/bin/python
# -*- coding: utf-8 -*-
import os
import glob
import shutil
from concurrent import futures
import datetime as dt
import threading
# 定义处理单个文件的函数
def process_file(file_path, to_folder_path):
'''脚本说明：
根据文件夹下的文件名，来新建文件夹，
并把所有该天的文件移动到新建的文件夹
'''
# 获取文件名
name = os.path.basename(file_path) # MYD04_L2.A2021001.0230.061.2021001233054.hdf
# 获取文件时间并转换为 %Y 年 %j 第几天
date_time = dt.datetime.strptime(name.split(".")[1], "A%Y%j")#A2021001
# 转换为自己想要的天的形式 这里依然利用 %Y 年 %j 第几天
day_name = dt.datetime.strftime(date_time, 'A%Y%j')#A2021001
# 定位到第几年第几日文件夹路径
folder_path = os.path.join(to_folder_path, str(date_time.year), day_name)
# 使用锁确保目录的原子性创建
with lock:
os.makedirs(folder_path, exist_ok=True)
# 之前文件的目录
old_file_path = file_path
# 新产生将要移动到的目录
new_file_path = os.path.join(folder_path, name)
# 移动文件（源文件路径，目标路径）
shutil.move(old_file_path, new_file_path)
if __name__ == '__main__':
# 需要整理的文件夹目录
folder_path = r"Y:\MODIS\AQUA\MYD04_10k_HDF"
# 需要整理到的文件夹目录
to_folder_path = r"Y:\MODIS\AQUA\MYD04_10k_HDF"
# 需要整理的hdf文件
files_list = glob.glob(os.path.join(folder_path, "*.hdf"))
# 使用锁确保线程安全
lock = threading.Lock()
# 使用ThreadPoolExecutor进行多线程处理
with futures.ThreadPoolExecutor() as executor:
# 并行处理每个文件
executor.map(process_file, files_list, [to_folder_path] * len(files_list))
print("所有文件处理完成")
文件夹整理之后，则为 Y:\MODIS\AQUA\MYD04_10k_HDF\2021\A2021001
二、时空匹配
1.对HDF文件进行匹配
#!/usr/bin/python
# -*- coding: utf-8 -*-
import os
import pandas as pd
import numpy as np
import datetime as dt
import glob
from tqdm import tqdm
from pyhdf.SD import SD, SDC
def read_site_csv(site_csv_path):
'''
@description: 读取站点数据，并根据指定的时间范围进行筛选
@param: 输入站点csv路径
@return: 返回读取站点数据dataframe
'''
#读取csv
aeronet_df=pd.read_csv(site_csv_path)
# 可更改站点起始终止范围时间
start_time = dt.datetime.strptime("2000-01-01 00:00", '%Y-%m-%d %H:%M')
end_time = dt.datetime.strptime("2021-01-01 00:00", '%Y-%m-%d %H:%M')
# 注意：这个函数假设 CSV 文件中有一个名为 "Date" 的列，且该列的日期格式为 'YYYY/MM/DD HH:MM'。如果实际数据的格式不同，你需要修改日期解析的代码。
aeronet_df["Date"]=aeronet_df["Date"].map(lambda x :dt.datetime.strptime(x, '%Y/%m/%d %H:%M'))
# 筛选在指定时间段内的数据
filtered_aeronet_df = aeronet_df[(aeronet_df['Date'] >= start_time) & (aeronet_df['Date'] <= end_time)]
return filtered_aeronet_df
def convert_Lat_Lon_to_row_col(lat, lon, lat_data, lon_data):
'''
@description: 将经纬度转换为 HDF 文件中的行列号
@param: lat, lon: 经纬度 lat_data, lon_data: HDF 文件中的经纬度
@return: row, col: HDF 文件中的行列号 min_distance: 最小距离(筛选作为判断依据)
'''
delta_lat = lat_data - lat
delta_lon = lon_data - lon
# 计算距离
distance = np.sqrt(delta_lat ** 2 + delta_lon ** 2)
# 获取最小距离的行列号
row, col = np.unravel_index(np.argmin(distance), distance.shape)
# 获取最小距离
min_distance = np.min(distance)
return row, col, min_distance
def hdf_geo_con(hdf_file_path,lon,lat):
'''
@description: 空间匹配判断经纬度是否在文件范围内
@param: hdf_file_path: HDF 文件路径 lon, lat: 经纬度
@return: True or False 范围是否在文件内
'''
try:
# 尝试打开 HDF 文件
hdf_file = SD(hdf_file_path, SDC.READ)
except Exception as e:
print(f"{hdf_file_path} 打开文件失败: {e}")
else:
#读取经纬度
datasets = hdf_file.datasets()
latarray = hdf_file.select('Latitude').get()
lonarray = hdf_file.select('Longitude').get()
# 获取经纬度范围
lon_min, lon_max = lonarray.min(),lonarray.max()
lat_min, lat_max = latarray.min(),latarray.max()
# 关闭 HDF 文件
hdf_file.end()
if (lon_min<=lon<=lon_max)&(lat_min<=lat<=lat_max):
return True
else:
return False
def read_hdf_data(hdf_file_path):
'''
@description: 读取HDF文件中的经纬度和AOD数据
@param: hdf_file_path HDF文件路径
@return: lon_data, lat_data, aod_data 经纬度和AOD数据
'''
try:
# 尝试打开 HDF 文件
hdf_file = SD(hdf_file_path, SDC.READ)
datasets = hdf_file.datasets()
lat_data = hdf_file.select('Latitude').get()
lon_data = hdf_file.select('Longitude').get()
aod_data = hdf_file.select('Deep_Blue_Aerosol_Optical_Depth_550_Land').get()
# 关闭 HDF 文件
hdf_file.end()
return lon_data, lat_data, aod_data
except HDF4Error as e:
print(f"无法打开 HDF 文件 {hdf_file_path}: {e}")
return None, None, None
def match_data(site_df,data_floder,data_suffix):
'''
@description: 站点数据与hdf数据时空匹配
@param: site_df 站点dataframe data_floder hdf数据文件夹 data_suffix hdf数据后缀
@return: df 匹配后的站点dataframe
'''
#复制站点数据,创建新的字段
df = site_df.copy()
df["aod550_11_pixel"] = None
df["aod550_33_pixel"] = None
df["min_distance"] = None
df["valid_pixels"] = None
#逐行匹配
for index,row in tqdm(df.iterrows(), total=df.shape[0]):
#对于每一行，通过列名访问对应的元素
#日期时间转换
date=dt.datetime.strptime(row["Date"], '%Y-%m-%d %H:%M:%S')
# date=dt.datetime.strptime(row["Date"], '%Y/%m/%d %H:%M') #2011/11/2 5:44
#获取经纬度
lon,lat = row["Site_Longitude(Degrees)"], row["Site_Latitude(Degrees)"]
#±30分钟匹配
date_start=date - dt.timedelta(minutes=30)
date_end=date + dt.timedelta(minutes=30)
#找到时间（天）对应的文件夹
#转换为 2020001（%Y年 %j一年当中的第几天）
date_time=dt.datetime.strftime(date,"A%Y%j")#A2006007
#找到数据当天文件夹
days_floder=os.path.join(data_floder,str(date.year),date_time) #../2006/A2006007
#定位到当天文件夹
datas=glob.glob(os.path.join(days_floder,data_suffix))
#对每个文件夹的文件进行循环匹配
for data in datas:
# 对路径名解析
data_name=os.path.basename(data) #MYD04_L2.A2020001.0355.061.2020002235451.hdf
pass_time=dt.datetime.strptime(data_name.split(".")[1]+data_name.split(".")[2],"A%Y%j%H%M")#A2020001.0355
# 判断当前时间是否在范围时间内
if date_start<=pass_time<=date_end:
#判断经纬度是否在文件范围内
if hdf_geo_con(data,lon,lat):
#读取经纬度和想要的数据
lon_array,lat_array,aod_array = read_hdf_data(data)
#未读取到数据退出
if aod_array is None:
continue
#找到经纬度对应的行列号
row, col, min_distance= convert_Lat_Lon_to_row_col(lat, lon, lat_array, lon_array)
#1*1像元单点匹配
aod_values_1_1=aod_array[row,col]
#筛除小于0的无效值
aod_values_1_1=np.where(aod_values_1_1<=0,np.nan,aod_values_1_1)
aod_values_1_1=np.nanmean(aod_values_1_1)/1000
# 匹配到数值
df.loc[index, 'aod550_11_pixel'] = aod_values_1_1
df.loc[index, 'min_distance'] = min_distance
#取3*3像元匹配
aod_values_3_3 = aod_array[row-1:row+2, col-1:col+2]
#筛除小于0的无效值
aod_values_3_3=np.where(aod_values_3_3<=0,np.nan,aod_values_3_3)
# 计算数组中的有效像元数量
valid_pixels = np.count_nonzero(~np.isnan(aod_values_3_3))
df.loc[index, 'valid_pixels']=valid_pixels
#3*3像元中 有效像元大于3个
if valid_pixels>=3:
#计算平均值
mean_aod = np.nanmean(aod_values_3_3)
#计算标准差
std_aod = np.nanstd(aod_values_3_3)
#符合正态分布筛除异常值
filtered_aod_values = [aod for aod in aod_values_3_3.flatten() if mean_aod - 3 * std_aod <= aod <= mean_aod + 3 * std_aod]
#筛除后有效像元大于3个
if len(filtered_aod_values) >= 3:
average_aod = np.nanmean(filtered_aod_values)/1000
# 使用 .loc[row_indexer, col_indexer] 设置数值
df.loc[index, 'aod550_33_pixel'] = average_aod
# print(f"{index} 匹配成功!")
return df
if __name__ == "__main__":
#站点csv 文件
site_csv_path=r"..\AOD20_Merger_CN_Polyfit2_AOD_2000_2020.csv"
#筛选需要读取站点数据
site_df=read_site_csv(site_csv_path)
#hdf文件夹
data_floder=r"..\C6.1\AQUA"
#hdf文件后缀
data_suffix="*.hdf"
#时空匹配
match_df=match_data(site_df,data_floder,data_suffix)
#筛选匹配到的数据
match_df=match_df[match_df['valid_pixels'] > 0]
#保存匹配到的数据
match_df.to_csv(r"..\AQUA_2000_2020_match.csv",index=False)
2.筛选结果
运行时使用了tqdm库显示进度条
匹配结果csv:
对匹配结果按需要筛除为匹配到的值，
1.筛除匹配到的空值
aod550_11_pixel 为空的值
2.筛除距离过大的值
半径25km ≈ 0.25°，min_distance<=0.25
3.筛除有效像元个数
筛除有效像元个数<3的值
结尾
文章主要利用AERONET 站点观测数据时空对比验证AOD数据产品，为验证提供一个思路，修改后续还可以更验证*.nc等数据，仅供参考。
希望以上内容能给大家带来帮助！
2302
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
