方案一：使用Arcpy處理
一、使用ArcMap處理
import arcpy
arcpy.env.overwriteOutput = True
inputnc = "E:/01UrbanWaterSourceSecurity/01Data/Runoff1992-2014/GRUN_v1_GSWP3_WGS84_05_1902_2014.nc"
# 提取出全部波段
arcpy.MakeNetCDFRasterLayer_md(inputnc,"Runoff","X","Y","Runoff_Layer","time","#","BY_VALUE") # 變量名稱等查看數據介紹，一般都會有或者參考方案二查看
# 提取單波段
i = 1 # 一般為從1開始
for yr in range(1902,2015): # 年份
for mon in range(1,13): # 月份
nctif = arcpy.MakeRasterLayer_management("Runoff_Layer","oneband","#","-180 -90 180 90",str(i)) # 提取單波段
name = "E:/01UrbanWaterSourceSecurity/01Data/Runoff1992-2014/tif/" + str(yr) + str(mon).zfill(2)+".tif"
arcpy.CopyRaster_management(nctif, name) # 單波段結果轉tif
print(name)
i = i + 1
方案二：使用python的netCDF4批量處理NC格式文件
一、使用ArcMap提取出第一期數據
1.使用工具箱中的“Make NetCDF Raster Layer”工具，提取出一個數據
可以發現該數據有正確的像元大小、坐標系等
2.導出該數據作為標准數據
二、使用ArcMap結合netCDF4
1. 查看數據屬性
from netCDF4 import Dataset,num2date
infile = "../01Data/Runoff1992-2014/GRUN_v1_GSWP3_WGS84_05_1902_2014.nc"
data_set = Dataset(infile) # 讀取nc文件信息
print(data_set)
輸出為
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
title: GRUN
version: GRUN 1.0
meteorological_forcing: GSWP3
temporal_resolution: monthly
spatial_resolution: 0.5x0.5
crs: WGS84
proj4: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
EPSG: 4326
references: Ghiggi et al.,2019. GRUN: An observation-based global gridded runoff dataset from 1902 to 2014. ESSD, doi: https://doi.org/10.5194/essd-2019-32
authors: Gionata Ghiggi; Lukas Gudmundsson
contacts: gionata.ghiggi@gmail.com; lukas.gudmundsson@env.ethz.ch
institution: Land-Climate Dynamics, Institute for Atmospheric and Climate Science, ETH Zürich
institution_id: IAC ETHZ
dimensions(sizes): X(720), Y(360), time(1356)
variables(dimensions): float64 X(X), float64 Y(Y), float64 time(time), float32 Runoff(time, Y, X)
groups:
可以看到variables變量X、Y為經緯度，time為時間，Runoff為需要的結果
2.批量導出結果
from osgeo import gdal
from netCDF4 import Dataset,num2date
import numpy as np
def WriteTiff(im_data,inputdir, path):
raster = gdal.Open(inputdir)
im_width = raster.RasterXSize #柵格矩陣的列數
im_height = raster.RasterYSize #柵格矩陣的行數
im_bands = raster.RasterCount #波段數
im_geotrans = raster.GetGeoTransform()#獲取仿射矩陣信息
im_proj = raster.GetProjection()#獲取投影信息
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
elif len(im_data.shape) == 2:
im_data = np.array([im_data])
else:
im_bands, (im_height, im_width) = 1, im_data.shape
# 創建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
if (dataset != None):
dataset.SetGeoTransform(im_geotrans) # 寫入仿射變換參數
dataset.SetProjection(im_proj) # 寫入投影
for i in range(im_bands):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset
infile = "../01Data/Runoff1992-2014/GRUN_v1_GSWP3_WGS84_05_1902_2014.nc"
data_set = Dataset(infile) # 讀取nc文件信息
time = data_set.variables["time"][:] # 獲取時間一列
units = data_set.variables["time"].units # 獲取第一期時間
#讀取樣本tif文件的地理信息
intif = "../03ProcessData/runoff_example.tif"
for i in range(0,len(time)):
yr = num2date(time[i],units).year # 提取年份
mon = num2date(time[i],units).month # 提取月份
value_data = data_set.variables['Runoff'][i]
# 將缺失值改為0
data = value_data.data
mask = value_data.mask
data[np.where(mask == True)] = 0
outputname = "../01Data/Runoff1992-2014/tif/" + str(yr) + str(mon).zfill(2) + ".tif"
WriteTiff(data,intif , outputname)
print(outputname)
方案三：只使用代碼netCDF4
import os
from osgeo import gdal,osr
from netCDF4 import Dataset,num2date
import numpy as np
import math
import matplotlib.pyplot as plt
#將程序、nc文件和樣本tif文件放在同一個目錄下
infile = r"D:\Data\pxh\grace\originalData\grace\CSR_GRACE_RL06_Mascons_all-corrections_v01.nc"
#讀取nc文件
os.chdir(os.path.abspath('.'))
data_set = Dataset(infile) # 讀取nc文件信息
time = data_set.variables["time"][:] # 獲取時間一列
units = data_set.variables["time"].units
latitude = data_set.variables["lat"][:] # 獲取緯度列
longitude = data_set.variables["lon"][:] # 獲取經度列
ind = next(x[0] for x in enumerate(longitude) if x[1] >= 180) # 獲取大於180的那一列的位置，后面用於從這個位置把后面的值移到前面解決東西半球反了的問題
# 為了提取出列表中的第一個值(即lon=180時所對應的index)
infileMask = r"D:\Data\pxh\grace\originalData\grace\CSR_GRACE_RL06_Mascons_v01_LandMask.nc"
data_set_Mask = Dataset(infileMask) # 讀取nc文件信息
value_mask = data_set_Mask.variables['LO_val'][:]
driver = gdal.GetDriverByName('GTiff')
dtype = gdal.GDT_Float32
osrs = osr.SpatialReference()
osrs.SetWellKnownGeogCS('WGS84')
osrs.ExportToWkt()
geoprojection = osrs.ExportToWkt()
print("開始運行")
for y in range(0,len(time)):
outdir = infile[:-3] # 獲取文件名，去除.nc
#建立文件夾，將生成的tif放在此文件夾中
if not os.path.isdir(outdir):
os.mkdir(outdir)
yr = num2date(time[y],units).year # 提取年份
mon = num2date(time[y],units).month # 提取月份
#讀取value data，將所得數組上下翻轉，lon大於180的部分挪到數組前面，方便最后直接寫入tif文件
value_data = data_set.variables['lwe_thickness'][:][y] # 獲取該時間對應的數據
# value_data = value_data[0,:,:] # 處理4維數據時使用
# value_data[np.where(value_mask == 0)] = np.nan # 海洋部分變為nan
value_data = np.flip(value_data, 0) # 矩陣上下左右反轉
value_data= np.roll(value_data,-ind,axis=1) # 從-ind-1位置把后面的值按列移到前面
#輸出tif文件
outtif = outdir + "/" + str(yr) + str(mon).zfill(2) +".tif"# 新建的文件夾/年份+月份（補充為2位的，1變為01等）+.tif
geotransform = (-180.0, 360.0/value_data.shape[1], 0.0, 90.0, 0.0, -180/value_data.shape[0]) # 360/重采樣后的列數得到輸出圖像的像元大小
outraster = driver.Create(outtif, value_data.shape[1],value_data.shape[0],1, dtype)
outraster.SetProjection(geoprojection)
outraster.SetGeoTransform(geotransform)
outraster.GetRasterBand(1).WriteArray(value_data)
outraster = None
print(outtif)
根據需求使用注釋掉的代碼
！注意事項
1.使用時候請自行修改修改輸入輸出文件路徑與變量名稱
2.根據需要處理缺失值
3.由於位置問題可能需要旋轉數據，推薦使用方案1
4.如果NC文件過大，使用arcmap處理過慢，可以使用方案3
溫度和降水數據月合成年
from osgeo import gdal,osr
import numpy as np
def LoadData(filename):
file = gdal.Open(filename)
if file == None:
print(filename + " can't be opened!")
return
nb = file.RasterCount
L = []
for i in range(1, nb + 1):
band = file.GetRasterBand(i)
background = band.GetNoDataValue()
data = band.ReadAsArray()
data = data.astype(np.float32)
index = np.where(data == background)
data[index] = 0
L.append(data)
data = np.stack(L,0)
if nb == 1:
data = data[0,:,:]
return data
def WriteTiff(im_data, im_width, im_height, im_bands, im_geotrans, im_proj, path):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
elif len(im_data.shape) == 2:
im_data = np.array([im_data])
else:
im_bands, (im_height, im_width) = 1, im_data.shape
# 創建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
if (dataset != None):
dataset.SetGeoTransform(im_geotrans) # 寫入仿射變換參數
dataset.SetProjection(im_proj) # 寫入投影
for i in range(im_bands):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset
indir = "TEM/tif/tmp_Layer2020_12.tif"
arr = LoadData(indir)
days = [31,28,31,30,31,30,31,31,30,31,30,31]
intifs = []
for i in range(len(days)):
intifs.append(arr[i] * days[i])
arr_mean = sum(intifs) / sum(days)
raster = gdal.Open(indir)
dtype = gdal.GDT_Float32
driver = gdal.GetDriverByName('GTiff')
im_width = raster.RasterXSize #柵格矩陣的列數
im_height = raster.RasterYSize #柵格矩陣的行數
im_bands = 1
osrs = osr.SpatialReference()
osrs.SetWellKnownGeogCS('WGS84')
osrs.ExportToWkt()
im_proj = osrs.ExportToWkt()
im_geotrans = raster.GetGeoTransform()#獲取仿射矩陣信息
ResultPath = "TEM/tif/tmp_Layer2020_mean.tif"
WriteTiff(arr_mean, im_width, im_height, im_bands, im_geotrans, im_proj, ResultPath)
#--------------------------------------------------------------------------
indir = "PRE/tif/pre_Layer1.tif"
arr = LoadData(indir)
intifs = []
for i in range(12):
intifs.append(arr[i])
arr_sum = sum(intifs)
raster = gdal.Open(indir)
dtype = gdal.GDT_Float32
driver = gdal.GetDriverByName('GTiff')
im_width = raster.RasterXSize #柵格矩陣的列數
im_height = raster.RasterYSize #柵格矩陣的行數
im_bands = 1
im_geotrans = raster.GetGeoTransform()#獲取仿射矩陣信息
osrs = osr.SpatialReference()
osrs.SetWellKnownGeogCS('WGS84')
osrs.ExportToWkt()
im_proj = osrs.ExportToWkt()
ResultPath = "PRE/tif/pre_Layer2020_sum1.tif"
WriteTiff(arr_sum, im_width, im_height, im_bands, im_geotrans, im_proj, ResultPath)