netCDF4的簡單操作

netcdf4是用處理netcdf檔案的函式庫

netcdf的全名是Network Common Data Form

副檔名縮寫為.nc

常用來儲存大氣、氣象、地理……等科學資料

netcdf的結構

netcdf可以分為4層架構:

1. group
2. Dimensions
3. Variables
4. Attributes

group

資料夾，用來分類資料

Dimensions

定義數據的時空間資訊，大小可固定可無限

Variables

實際數據，必須綁定一個Dimensions，可以附帶Attributes

Attributes

描述group、Dimensions、Variables的資訊，如:單位、文字描述

netCDF4的crud

create

創建nc檔案

rootgrp = Dataset('test.nc', 'w', format='NETCDF4')
print(rootgrp.file_format)
rootgrp.close()

創建group

rootgrp = Dataset('test.nc', 'a')
fcstgrp = rootgrp.createGroup('forecasts')
analgrp = rootgrp.createGroup('analyses')
fcstgrp1 = rootgrp.createGroup('/forecasts/model1')
fcstgrp2 = rootgrp.createGroup('/forecasts/model2')

print(rootgrp.groups)

創建Dimensions

level_dim = rootgrp.createDimension('level', None)
time_dim = rootgrp.createDimension('time', None)
lat_dim = rootgrp.createDimension('lat', 73)
lon_dim = rootgrp.createDimension('lon', 144)

第二個參數就是決定Dimensions的大小，假設輸入特定數字，代表Dimensions具有大小，反之，Dimensions的大小沒有限制，可以動態增加。

創建Variables

times = rootgrp.createVariable('time','f8',('time',))
levels = rootgrp.createVariable('level','i4',('level',))
latitudes = rootgrp.createVariable('lat','f4',('lat',))
longitudes = rootgrp.createVariable('lon','f4',('lon',))
#temp包含四個維度，並進行壓縮處理，只保留小數點後三位的精度
temp = rootgrp.createVariable('temp','f4',('time','level','lat','lon',),least_significant_digit=3)
#在group創建變數temp
temp = rootgrp.createVariable('/forecasts/model1/temp','f4',('time','level','lat','lon',))

創建Attributes

import time
rootgrp.description = 'bogus example script'
rootgrp.history = 'Created ' + time.ctime(time.time())
rootgrp.source = 'netCDF4 python module tutorial'
latitudes.units = 'degrees north'
longitudes.units = 'degrees east'
levels.units = 'hPa'
temp.units = 'K'
times.units = 'hours since 0001-01-01 00:00:00.0'
calendar: Literal['gregorian'] = 'gregorian'
times.calendar = calendar

for name in rootgrp.ncattrs():
    print('Global attr', name, '=', getattr(rootgrp,name))

print(rootgrp)

print(rootgrp.__dict__)

print(rootgrp.variables)

對Variables賦值

lats = np.arange(-90, 91, 2.5)
lons = np.arange(-180, 180, 2.5)
latitudes[:] = lats
longitudes[:] = lons 

# append along two unlimited dimensions by assigning to slice.
nlats = len(rootgrp.dimensions['lat'])
nlons = len(rootgrp.dimensions['lon'])
print('temp shape before adding data = ',temp.shape)
# random number generator.
temp[0:5,0:10,:,:] = uniform(size=(5,10,nlats,nlons))
print('temp shape after adding data = ',temp.shape)
# levels have grown, but no values yet assigned.
print('levels shape after adding pressure data = ',levels.shape)

dates = [datetime(2001,3,1)+n*timedelta(hours=12) for n in range(temp.shape[0])]
times[:] = date2num(dates,units=times.units,calendar=times.calendar)
print("time values (in units {}):\n{}".format(times.units, times[:]))
dates_array = num2date(times[:],units=times.units,calendar=times.calendar)
print("dates corresponding to time values:\n{}".format(dates_array))

對複合類型(compound type)賦值

# more complex compound type example.
nc = Dataset('compound_example.nc','w') # create a new dataset.
# create an unlimited  dimension call 'station'
nc.createDimension('station',None)
# define a compound data type (can contain arrays, or nested compound types).
winddtype = np.dtype([('speed','f4'),('direction','i4')])
statdtype = np.dtype([('latitude', 'f4'), ('longitude', 'f4'),
                      ('surface_wind',winddtype),
                      ('temp_sounding','f4',10),('press_sounding','i4',10),
                      ('location_name','S12')])
# use this data type definitions to create a compound data types
# called using the createCompoundType Dataset method.
# create a compound type for vector wind which will be nested inside
# the station data type. This must be done first!
wind_data_t = nc.createCompoundType(winddtype,'wind_data')
# now that wind_data_t is defined, create the station data type.
station_data_t = nc.createCompoundType(statdtype,'station_data')
# create nested compound data types to hold the units variable attribute.
winddtype_units = np.dtype([('speed','S12'),('direction','S12')])
statdtype_units = np.dtype([('latitude', 'S12'), ('longitude', 'S12'),
                            ('surface_wind',winddtype_units),
                            ('temp_sounding','S12'),
                            ('location_name','S12'),
                            ('press_sounding','S12')])
# create the wind_data_units type first, since it will nested inside
# the station_data_units data type.
wind_data_units_t = nc.createCompoundType(winddtype_units,'wind_data_units')
station_data_units_t =\
nc.createCompoundType(statdtype_units,'station_data_units')
# create a variable of of type 'station_data_t'
statdat = nc.createVariable('station_obs', station_data_t, ('station',))
# create a numpy structured array, assign data to it.
data = np.empty(1,statdtype)
data['latitude'] = 40.
data['longitude'] = -105.
data['surface_wind']['speed'] = 12.5
data['surface_wind']['direction'] = 270
data['temp_sounding'] = (280.3,272.,270.,269.,266.,258.,254.1,250.,245.5,240.)
data['press_sounding'] = range(800,300,-50)
data['location_name'] = 'Boulder, CO'
# assign structured array to variable slice.
statdat[0] = data
# or just assign a tuple of values to variable slice
# (will automatically be converted to a structured array).
statdat[1] = np.array((40.78,-73.99,(-12.5,90),
             (290.2,282.5,279.,277.9,276.,266.,264.1,260.,255.5,243.),
             range(900,400,-50),'New York, NY'),data.dtype)
print(nc.cmptypes)
windunits = np.empty(1,winddtype_units)
stationobs_units = np.empty(1,statdtype_units)
windunits['speed'] = 'm/s'
windunits['direction'] = 'degrees'
stationobs_units['latitude'] = 'degrees N'
stationobs_units['longitude'] = 'degrees W'
stationobs_units['surface_wind'] = windunits
stationobs_units['location_name'] = 'None'
stationobs_units['temp_sounding'] = 'Kelvin'
stationobs_units['press_sounding'] = 'hPa'
print(stationobs_units.dtype)
statdat.units = stationobs_units
# close and reopen the file.
nc.close()
nc = Dataset('compound_example.nc')
print(nc)
statdat = nc.variables['station_obs']
print(statdat)
# print out data in variable.
print('data in a variable of compound type:')
print(statdat[:])
nc.close()

上面的程式是要記錄氣象站的資料

dtype是Variables的屬性之一，紀錄Variables的資料型態

上面times = rootgrp.createVariable(‘time’,’f8',(‘time’,))的'f8'就是dtype

之所以將dtype獨立出來，是因為這裡遇到複合類型:這是指一個變數由多個字段組成，每個字段有自己的資料型態

np.dtype就是定義各個字段的資料型態

其他部分與一般賦值相同，建立Variables or Attributes，之後輸入數值

read

查看特定Variables

print(f.variables.keys()) # get all variable names
temp = f.variables['temperature']  # temperature variable
print(temp)

查看Dimensions

for d in f.dimensions.items():
    print(d)

查看特定Variables的Dimensions與大小

temp.dimensions
temp.shape

查找特定經緯度的資訊

lat, lon = f.variables['Latitude'], f.variables['Longitude']
print(lat)

# extract lat/lon values (in degrees) to numpy arrays
latvals = lat[:]; lonvals = lon[:] 
# a function to find the index of the point closest pt
# (in squared distance) to give lat/lon value.
def getclosest_ij(lats,lons,latpt,lonpt):
    # find squared distance of every point on grid
    dist_sq = (lats-latpt)**2 + (lons-lonpt)**2  
    # 1D index of minimum dist_sq element
    minindex_flattened = dist_sq.argmin()    
    # Get 2D index for latvals and lonvals arrays from 1D index
    return np.unravel_index(minindex_flattened, lats.shape)
iy_min, ix_min = getclosest_ij(latvals, lonvals, 50., -140)

缺失值

soilmvar = gfs.variables['Volumetric_Soil_Moisture_Content_depth_below_surface_layer']
# flip the data in latitude so North Hemisphere is up on the plot
soilm = soilmvar[0,0,::-1,:] 
print('shape=%s, type=%s, missing_value=%s' % \
      (soilm.shape, type(soilm), soilmvar.missing_value))
import matplotlib.pyplot as plt
%matplotlib inline
cs = plt.contourf(soilm)

netcdf中，Variables有一個屬性表示缺失值，NetCDF4讀到缺失值時，會將變數轉為MaskedArray

轉成MaskedArray後，NetCDF4會自動忽略缺失值

當我們將資料可視化，白色部分就是缺失值

update

使用前面複合類型的例子，假如我們想把第一筆資料的地名改成台北，

我們可以透過以下操作完成任務:

nc = Dataset('compound_example.nc','r+')
statdat = nc.variables['station_obs']
data = statdat[0]
data['location_name'] = 'taipei'
statdat[0] = data
print(statdat[0])

nc.close()

delete

文件內沒翻到delete，如果有熟悉這塊的大佬，歡迎在下面補充。

參考資料:

netcdf4-python/examples/tutorial.py at master · Unidata/netcdf4-python

netcdf4-python/examples/reading_netCDF.ipynb at master · Unidata/netcdf4-python