数据为全球近20年的气温及总降水量NetCDF(nc)影像,处理为tif格式的影响,并裁剪为中国区域。
import arcpy
from arcpy import env
from arcpy.sa import *
#t2m对应温度波段、tp对应降水波段
variable = "tp" #对应的波段(不知道是不是叫波段,因为一张影像里面既有温度也有降水)
XDimension = "longitude"
YDimension = "latitude"
#outRasterLayer = "t2m_layer"
bandDimmension = ""
dimensionValues = ""
valueSelectionMethod = ""
cellRegistration = ""
for year in range(2001, 2021): # 年份
env.workspace = r"Z:\VM\NDVI2\tem&pre" + '\\' + str(year) #输入的数据路径
outpath = r"Z:\VM\NDVI2\tem&pre_mask\pre" + '\\' + str(year) # 输出的数据路径
#rasters = arcpy.ListRasters("*", All) #无效,arcpy获取不了nc文件
rasters = [] #用for循环来获取所有nc文件
for i in range(1,13):
if i < 10:
rasters.append(r"Z:\VM\NDVI2\tem&pre" + '\\' + str(year) + '\\' + str(year) + "-0" + str(i) + "tem&pre.nc") #按需修改
else:
rasters.append(r"Z:\VM\NDVI2\tem&pre" + '\\' + str(year) + '\\' + str(year) + "-" + str(i) + "tem&pre.nc") #同上
#value是下面Minus函数用的,为了将华氏温度转为摄氏度(-273.15)
#因为Minus的第二个参数放数值一直报错,所以我就生成了一张所有像元(像素)值都为273.15的tif影像
#生成方式:见代码下方
#value = r"Z:\VM\NDVI2\27315\27315.tif"
for raster in rasters: # 遍历影像
out = outpath+'\\'+(raster.split("-", 1)[1]).split("tem&",1)[0] +((raster.split("-", 1)[1]).split("tem&",1)[1]).split(".",1)[0] + ".tif" #输出位置及名称
arcpy.MakeNetCDFRasterLayer_md(raster, variable, XDimension, YDimension,
out, bandDimmension, dimensionValues,
valueSelectionMethod, cellRegistration)
#outMinus = Minus(out,value) #减函数
outExtractByMask = ExtractByMask(out, r"Z:\VM\NDVI2\map\bou2_4p.shp") # 掩膜、去掉影像的黑边,更美观
outExtractByMask.save(out)
生成方式:
arcgis先导入一个nc文件,「多维工具-创建NetCDF栅格图册」,变量改成需要的,其它默认
「spatial analyst-条件分析-条件函数」,不管true or false全部都为273.15