用于处理 NASA AOD 数据/可视化的脚本

nasa_aod.py

import time

import h5py
import matplotlib.pyplot as plt
import numpy as np


def timer(func):
    """给函数计时的装饰器"""

    def wrapped(*args, **kwargs):
        start_time = time.time()
        ret = func(*args, **kwargs)
        end_time = time.time()
        print(f"{func.__name__}: 用时 {end_time - start_time:.3f}s.")
        return ret

    return wrapped


@timer
def calibrate(aod_data, offset, scale, fill_val, valid_min, valid_max):
    """用于对数据进行校正，返回校正后的数据，以及有效数据的 mask（mask 为 True 的为有效数据）"""
    invalid_mask = (
        (aod_data < valid_min) | (aod_data > valid_max) | (aod_data == fill_val)
    )
    aod_data = (aod_data - offset) * scale
    return aod_data, ~invalid_mask


@timer
def visualize(aod_data, masks):
    """用于对校正后的数据进行可视化"""
    coo_vec = np.arange(1200)
    xx, yy = np.meshgrid(coo_vec, coo_vec, indexing="ij")

    for orbit in range(aod_data.shape[0]):
        channel = aod_data[orbit]
        mask = masks[orbit]
        plt.scatter(
            xx[mask].ravel(),
            yy[mask].ravel(),
            s=0.05,
            c=channel[mask].ravel(),
            cmap="gist_rainbow_r",
        )
    plt.colorbar()
    plt.axis("equal")


if __name__ == "__main__":

    # 打开文件，拿到 AOD 数据
    h5file = h5py.File("MCD19A2.A2022240.h15v05.006.2022242033236.h5", "r")
    Optical_Depth_047 = h5file["/grid1km/Data Fields/Optical_Depth_047"]

    # 从数据的属性中提取校正所需的参数
    add_offset = Optical_Depth_047.attrs["add_offset"]
    scale_factor = Optical_Depth_047.attrs["scale_factor"]
    fill_value = Optical_Depth_047.attrs["_FillValue"]
    valid_min, valid_max = Optical_Depth_047.attrs["valid_range"]
    print(
        f"校正参数：add_offset: {add_offset}, scale_factor: {scale_factor}, fill_value: {fill_value},",
        f"valid_range: [{valid_min}, {valid_max}]",
    )

    # 先从 int 转成 double，便于后面校正计算
    aod_data: np.ndarray = Optical_Depth_047[...].astype(float)

    # 对数据进行校正
    aod_data, masks = calibrate(
        aod_data, add_offset, scale_factor, fill_value, valid_min, valid_max
    )

    # print("校正后：")
    # print(aod_data)

    visualize(aod_data, masks)

    # 显示可视化结果窗口
    plt.show()