快捷導(dǎo)航

python中的netCDF4批量處理NC文件的操作方法

更新時(shí)間：2022年03月21日 15:01:38 作者：skypanxh

這篇文章主要介紹了python的netCDF4批量處理NC格式文件的操作方法,使用python批量提取所有數(shù)據(jù)，查看數(shù)據(jù)屬性，本文通過(guò)實(shí)例代碼給大家介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或工作具有一定的參考借鑒價(jià)值，需要的朋友可以參考下

一、使用ArcMap提取出第一期數(shù)據(jù)

1.使用工具箱中的“Make NetCDF Raster Layer”工具，提取出一個(gè)數(shù)據(jù)

可以發(fā)現(xiàn)該數(shù)據(jù)有正確的像元大小、坐標(biāo)系等

2.導(dǎo)出該數(shù)據(jù)作為標(biāo)準(zhǔn)數(shù)據(jù)

二、使用python批量提取所有數(shù)據(jù)

1. 查看數(shù)據(jù)屬性

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為經(jīng)緯度，time為時(shí)間，Runoff為需要的結(jié)果

2.批量導(dǎo)出結(jié)果

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 #柵格矩陣的列數(shù)
    im_height = raster.RasterYSize #柵格矩陣的行數(shù)
    im_bands = raster.RasterCount #波段數(shù)
    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])
        im_bands, (im_height, im_width) = 1, im_data.shape
        # 創(chuàng)建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
    if (dataset != None):
        dataset.SetGeoTransform(im_geotrans)  # 寫(xiě)入仿射變換參數(shù)
        dataset.SetProjection(im_proj)  # 寫(xiě)入投影
    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"][:]  # 獲取時(shí)間一列
units = data_set.variables["time"].units # 獲取第一期時(shí)間
#讀取樣本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)