利用Python实现对葵花8号L1级反射率数据进行云格点剔除的操作实例

一、内容介绍

本文利用python下的gdal库，对葵花8号影像L1级反射率数据，基于L2级的云检测产品数据
（CLP）下的云类型对有云格点数据进行筛选，剔除操作，得到剔除云格点的反射率数据。
本文的反射率数据和云类型数据均为tif格式，需要将原始数据nc格式转为tif格式。

二、代码分块讲解

1、导入所需的库

import os
import numpy as np
from osgeo import gdal
import glob

2、输入与输出路径设置

AOD_file_path = r"F:\PM2.5反演代码\test_cloud\AOD"
cloud_file_path = r"F:\PM2.5反演代码\test_cloud\cloud_type"
out_put = r"F:\PM2.5反演代码\test_cloud\AOD_clear"

3、读取输入路径下的反射率数据和云类型数据，数据格式均为tif

AOD_file_list =  glob.glob(AOD_file_path + '/*.tif')
cloud_file_list = glob.glob(cloud_file_path + '/*.tif')

4、依据命名匹配对应反射率下的云类型数据

for AOD_file in AOD_file_list:   #循环读取
    AOD_name = AOD_file.split("\\")[-1].split("_")[0]

    for cd in cloud_file_list:
        cd_name = cd.split("\\")[-1].split("_")[0]
        if AOD_name in cd_name:
            AOD =  AOD_file
            cloud =  cd

5、读取反射率数据和云格点数据为矩阵格式

#AOT数据读取
AOD_raster = gdal.Open(AOD)
AOD_array = AOD_raster.ReadAsArray()
#cloud数据读取
cd_raster = gdal.Open(cloud)
cd_array = cd_raster.ReadAsArray()

6、对云格点数据进行剔除（葵花8号影像云格点数据clear = 0,其他为有云）

#质量筛选
AOD_result = np.where(cd_array != 0,np.nan,AOD_array)

7、读取影像参数

#影像参数读取
row = AOD_raster.RasterYSize
col = AOD_raster.RasterXSize
geotransform = AOD_raster.GetGeoTransform()
projection=AOD_raster.GetProjection()

8、保存结果

#保存结果
pos = out_put + os.sep + AOD_name + ".tif"
driver=gdal.GetDriverByName("Gtiff")
result = driver.Create(pos,row,col,1,gdal.GDT_Float32)
result.SetGeoTransform(geotransform)
result.SetProjection(projection)
result.GetRasterBand(1).WriteArray(AOD_result)
result.GetRasterBand(1).SetNoDataValue(eval("-32768"))
result.FlushCache()  # 将数据写入硬盘
del result  # 注意必须关闭tif文件

三、完整代码

import os
import numpy as np
from osgeo import gdal
import glob

AOD_file_path = r"F:\PM2.5反演代码\test_cloud\AOD"
cloud_file_path = r"F:\PM2.5反演代码\test_cloud\cloud_type"
out_put = r"F:\PM2.5反演代码\test_cloud\AOD_clear"

AOD_file_list =  glob.glob(AOD_file_path + '/*.tif')
cloud_file_list = glob.glob(cloud_file_path + '/*.tif')

for AOD_file in AOD_file_list:
    AOD_name = AOD_file.split("\\")[-1].split("_")[0]

    for cd in cloud_file_list:
        cd_name = cd.split("\\")[-1].split("_")[0]
        if AOD_name in cd_name:

            AOD =  AOD_file
            cloud =  cd

            #AOD数据读取
            AOD_raster = gdal.Open(AOD)
            AOD_array = AOD_raster.ReadAsArray()
            #cloud数据读取
            cd_raster = gdal.Open(cloud)
            cd_array = cd_raster.ReadAsArray()

            #影像参数读取
            row = AOD_raster.RasterYSize
            col = AOD_raster.RasterXSize
            geotransform = AOD_raster.GetGeoTransform()
            projection=AOD_raster.GetProjection()

            #质量筛选
            AOD_result = np.where(cd_array != 0,np.nan,AOD_array)

            #保存结果
            pos = out_put + os.sep + AOD_name + ".tif"
            driver=gdal.GetDriverByName("Gtiff")
            result = driver.Create(pos,row,col,1,gdal.GDT_Float32)
            result.SetGeoTransform(geotransform)
            result.SetProjection(projection)
            result.GetRasterBand(1).WriteArray(AOD_result)
            result.GetRasterBand(1).SetNoDataValue(eval("-32768"))
            result.FlushCache()  # 将数据写入硬盘
            del result  # 注意必须关闭tif文件

四、参考文章

https://www.cnblogs.com/fkxxgis/p/17147496.html