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

大美女

两年前，我曾发布过一篇名为《导出 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 代码：

1. import os
2. import arcpy
4. # 脚本工具核心代码
5. def script_tool(nc_folder, variable, x_dimension, y_dimension, band_dimension, tif_folder):
6.     # 基于工作空间为 tif_folder，后续输出文件都将存放在该文件夹中
7.     arcpy.env.workspace = tif_folder
8.     # 遍历文件夹中所有后缀名为 .nc 的文件，将其逐一转换为 tif 文件
9.     for file in os.listdir(nc_folder):
10.         if file[-3:] == ".nc":
11.             # 以原始文件名称为基础，得到不包含后缀名部分的 tif 文件名
12.             tif_file_name = file[:-3]
13.             # 判断 band_dimension 参数是否存在输入值，有值则在后续函数中指定 band_dimension 参数
14.             if band_dimension:
15.                 # 基于 NetCDF 文件创建栅格图层
16.                 arcpy.MakeNetCDFRasterLayer_md(f"{nc_folder}/{file}", variable, x_dimension, y_dimension, tif_file_name, band_dimension)
17.             else:
18.                 arcpy.MakeNetCDFRasterLayer_md(f"{nc_folder}/{file}", variable, x_dimension, y_dimension, tif_file_name)
19.             # 基于栅格图层创建 tif 文件
20.             arcpy.CopyRaster_management(tif_file_name, tif_file_name + ".tif")
21.             arcpy.AddMessage(tif_file_name + " " + "conversion successful!")
22.     return
24. # 读取输入参数，执行程序
25. if __name__ == "__main__":
26.     # 获取 ArcGIS 脚本工具参数
27.     nc_folder = arcpy.GetParameterAsText(0)
28.     variable = arcpy.GetParameterAsText(1)
29.     x_dimension = arcpy.GetParameterAsText(2)
30.     y_dimension = arcpy.GetParameterAsText(3)
31.     band_dimension = arcpy.GetParameterAsText(4)
32.     tif_folder = arcpy.GetParameterAsText(5)
33.     # 执行脚本工具核心代码
34.     script_tool(nc_folder, variable, x_dimension, y_dimension, band_dimension, tif_folder)

复制代码

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

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

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

2.2 遍历并导出栅格各波段

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

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

1. import os
2. import arcpy
4. def script_tool(input_folder, output_folder):
5.     # 遍历文件夹中所有 tif 文件
6.     for file in os.listdir(input_folder):
7.         if file[-4:] == ".tif":
8.             # 根据文件名称获取对应年份，此处为倒数第5位至倒数第8位，可根据实际情况调整
9.             year = file[-8:-4]
10.             # 创建循环，遍历各月份对应波段，将波段逐一导出为单独的 tif 文件并保存到输出文件夹
11.             for i in range(12):
12.                 input_band = f"{input_folder}/{file}/Band_{i+1}"
13.                 output_raster = f"{output_folder}/Pre_{year}_{i + 1}.tif"
14.                 arcpy.CopyRaster_management(input_band, output_raster)
15.                 arcpy.AddMessage(f"Pre_{year}_{i+1} conversion successful!")
16.     return
18. if __name__ == "__main__":
19.     # 获取 ArcGIS 脚本工具参数
20.     input_folder = arcpy.GetParameterAsText(0)
21.     output_folder = arcpy.GetParameterAsText(1)
22.     # 执行脚本工具核心代码
23.     script_tool(input_folder, output_folder)

复制代码

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

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

来源:https://www.cnblogs.com/qsgeo/p/18354723
免责声明：由于采集信息均来自互联网，如果侵犯了您的权益，请联系我们【E-Mail:cb@itdo.tech】 我们会及时删除侵权内容，谢谢合作！