NetCDF 文件批量转栅格并导出栅格各波段

2024-08-12 12:44 由姜颢睿发表于 #后端开发

两年前，我曾发布过一篇名为《导出 NetCDF 栅格图层的各个波段》的公众号推文，讲述了通过网络中的ArcGIS工具将单个 NetCDF 文件的各个波段分别导出为 tif 文件的方法。该工具提供了 arcpy 源代码，我们以该代码为基础，将其转换为 ArcGIS Pro 环境下的 Python 3 代码，并使程序可对多个文件进行批处理，本文将对该代码及其在 ArcGIS Pro 中的配置过程进行介绍。

1 数据来源及介绍

本文所使用的数据为国家青藏高原科学数据中心的中国1km分辨率逐月降水量数据集。对下载得到的 nc 文件通过 Panoply 软件进行读取可知，其X、Y及波段维度名称分别为lon、lat、time，nc文件中共12个波段，分别存储该年度1-12月每个月份的降水量。

2 脚本配置及相关代码

2.1 NetCDF 转换为 tif

在 ArcGIS Pro 地图工程的默认工具箱中新建脚本工具，命名为 NetCDF to TIF，并仿照 ArcGIS 的 创建 NetCDF 栅格图层 工具设置其参数。该工具用于读取文件夹中的所有 NetCDF 文件，根据其 X、Y、波段等维度将其转换为 TIF 文件并逐一输出至 TIF 文件夹中。由于单波段 NetCDF 文件不存在波段维度，故该参数设置为可选类型。

右键单击脚本工具，在弹出菜单中单击编辑，将编辑窗口中的默认代码替换为下方 Python 代码：

import os
import arcpy

# 脚本工具核心代码
def script_tool(nc_folder, variable, x_dimension, y_dimension, band_dimension, tif_folder):
    # 基于工作空间为 tif_folder，后续输出文件都将存放在该文件夹中
    arcpy.env.workspace = tif_folder
    # 遍历文件夹中所有后缀名为 .nc 的文件，将其逐一转换为 tif 文件
    for file in os.listdir(nc_folder):
        if file[-3:] == ".nc":
            # 以原始文件名称为基础，得到不包含后缀名部分的 tif 文件名
            tif_file_name = file[:-3]
            # 判断 band_dimension 参数是否存在输入值，有值则在后续函数中指定 band_dimension 参数
            if band_dimension:
                # 基于 NetCDF 文件创建栅格图层
                arcpy.MakeNetCDFRasterLayer_md(f"{nc_folder}/{file}", variable, x_dimension, y_dimension, tif_file_name, band_dimension)
            else:
                arcpy.MakeNetCDFRasterLayer_md(f"{nc_folder}/{file}", variable, x_dimension, y_dimension, tif_file_name)
            # 基于栅格图层创建 tif 文件
            arcpy.CopyRaster_management(tif_file_name, tif_file_name + ".tif")
            arcpy.AddMessage(tif_file_name + " " + "conversion successful!")
    return

# 读取输入参数，执行程序
if __name__ == "__main__":
    # 获取 ArcGIS 脚本工具参数
    nc_folder = arcpy.GetParameterAsText(0)
    variable = arcpy.GetParameterAsText(1)
    x_dimension = arcpy.GetParameterAsText(2)
    y_dimension = arcpy.GetParameterAsText(3)
    band_dimension = arcpy.GetParameterAsText(4)
    tif_folder = arcpy.GetParameterAsText(5)
    # 执行脚本工具核心代码
    script_tool(nc_folder, variable, x_dimension, y_dimension, band_dimension, tif_folder)

保存后可在右键菜单的工具属性中看到下方界面：

配置完成后，将 Panoply 读取到的参数输入脚本，如下图所示：

程序运行完成后输出的栅格数据如下图所示：

2.2 遍历并导出栅格各波段

新建脚本工具，命名为 Band Split，并设置参数。

将脚本工具的默认代码替换为下方 Python 代码：

import os
import arcpy

def script_tool(input_folder, output_folder):
    # 遍历文件夹中所有 tif 文件
    for file in os.listdir(input_folder):
        if file[-4:] == ".tif":
            # 根据文件名称获取对应年份，此处为倒数第5位至倒数第8位，可根据实际情况调整
            year = file[-8:-4]
            # 创建循环，遍历各月份对应波段，将波段逐一导出为单独的 tif 文件并保存到输出文件夹
            for i in range(12):
                input_band = f"{input_folder}/{file}/Band_{i+1}"
                output_raster = f"{output_folder}/Pre_{year}_{i + 1}.tif"
                arcpy.CopyRaster_management(input_band, output_raster)
                arcpy.AddMessage(f"Pre_{year}_{i+1} conversion successful!")
    return

if __name__ == "__main__":
    # 获取 ArcGIS 脚本工具参数
    input_folder = arcpy.GetParameterAsText(0)
    output_folder = arcpy.GetParameterAsText(1)
    # 执行脚本工具核心代码
    script_tool(input_folder, output_folder)

配置完成后，将上一步得到的 tif 文件夹作为输入文件夹参数输入脚本，如下图所示：

程序运行完成后输出的栅格数据如下图所示：