翼度科技»论坛 编程开发 python 查看内容

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

5

主题

5

帖子

15

积分

新手上路

Rank: 1

积分
15
两年前,我曾发布过一篇名为《导出 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
  3. # 脚本工具核心代码
  4. def script_tool(nc_folder, variable, x_dimension, y_dimension, band_dimension, tif_folder):
  5.     # 基于工作空间为 tif_folder,后续输出文件都将存放在该文件夹中
  6.     arcpy.env.workspace = tif_folder
  7.     # 遍历文件夹中所有后缀名为 .nc 的文件,将其逐一转换为 tif 文件
  8.     for file in os.listdir(nc_folder):
  9.         if file[-3:] == ".nc":
  10.             # 以原始文件名称为基础,得到不包含后缀名部分的 tif 文件名
  11.             tif_file_name = file[:-3]
  12.             # 判断 band_dimension 参数是否存在输入值,有值则在后续函数中指定 band_dimension 参数
  13.             if band_dimension:
  14.                 # 基于 NetCDF 文件创建栅格图层
  15.                 arcpy.MakeNetCDFRasterLayer_md(f"{nc_folder}/{file}", variable, x_dimension, y_dimension, tif_file_name, band_dimension)
  16.             else:
  17.                 arcpy.MakeNetCDFRasterLayer_md(f"{nc_folder}/{file}", variable, x_dimension, y_dimension, tif_file_name)
  18.             # 基于栅格图层创建 tif 文件
  19.             arcpy.CopyRaster_management(tif_file_name, tif_file_name + ".tif")
  20.             arcpy.AddMessage(tif_file_name + " " + "conversion successful!")
  21.     return
  22. # 读取输入参数,执行程序
  23. if __name__ == "__main__":
  24.     # 获取 ArcGIS 脚本工具参数
  25.     nc_folder = arcpy.GetParameterAsText(0)
  26.     variable = arcpy.GetParameterAsText(1)
  27.     x_dimension = arcpy.GetParameterAsText(2)
  28.     y_dimension = arcpy.GetParameterAsText(3)
  29.     band_dimension = arcpy.GetParameterAsText(4)
  30.     tif_folder = arcpy.GetParameterAsText(5)
  31.     # 执行脚本工具核心代码
  32.     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
  3. def script_tool(input_folder, output_folder):
  4.     # 遍历文件夹中所有 tif 文件
  5.     for file in os.listdir(input_folder):
  6.         if file[-4:] == ".tif":
  7.             # 根据文件名称获取对应年份,此处为倒数第5位至倒数第8位,可根据实际情况调整
  8.             year = file[-8:-4]
  9.             # 创建循环,遍历各月份对应波段,将波段逐一导出为单独的 tif 文件并保存到输出文件夹
  10.             for i in range(12):
  11.                 input_band = f"{input_folder}/{file}/Band_{i+1}"
  12.                 output_raster = f"{output_folder}/Pre_{year}_{i + 1}.tif"
  13.                 arcpy.CopyRaster_management(input_band, output_raster)
  14.                 arcpy.AddMessage(f"Pre_{year}_{i+1} conversion successful!")
  15.     return
  16. if __name__ == "__main__":
  17.     # 获取 ArcGIS 脚本工具参数
  18.     input_folder = arcpy.GetParameterAsText(0)
  19.     output_folder = arcpy.GetParameterAsText(1)
  20.     # 执行脚本工具核心代码
  21.     script_tool(input_folder, output_folder)
复制代码
配置完成后,将上一步得到的 tif 文件夹作为输入文件夹参数输入脚本,如下图所示:

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


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

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?立即注册

x
来自手机

举报 回复 使用道具