点击上方 
蓝字关注我们
气象领域很多数据都是以nc格式存储的，这类数据同时封装了数据和描述信息。Python、NCL、Fortran和Matlab等软件都可以对其进行读写操作，这里记录一下利用netCDF4和xarray两个库读写nc数据。
import netCDF4 as nc
filename = "FCST_d01_2018-07-27_12_00_00.nc"fn = nc.Dataset(filename,"r")
print(fn)'''root group (NETCDF3_CLASSIC data model, file format NETCDF3): creation_date: 2019年 11月 02日 星期六 21:59:04 CST NCL_Version: 6.4.0 system: Linux rlk.iap 3.10.0-957.12.1.el7.x86_64 #1 SMP Mon Apr 29 \ 14:59:59 UTC 2019 x86_64 x86_64 x86_64 GNU/Linux Conventions: None grib_source: FCST_d01_2018-07-27_12_00_00.grib title: NCL: convert-GRIB-to-netCDF dimensions(sizes): g1_lat_0(921), g1_lon_1(804) variables(dimensions): float32 TMP_GDS1_HTGL(g1_lat_0,g1_lon_1), float32 VIS_GDS1_SFC(g1_lat_0,g1_lon_1), \ float32 R_H_GDS1_HTGL(g1_lat_0,g1_lon_1), float32 A_PCP_GDS1_SFC_acc12h(g1_lat_0,g1_lon_1), \ float32 NCPCP_GDS1_SFC_acc12h(g1_lat_0,g1_lon_1), float32 ACPCP_GDS1_SFC_acc12h(g1_lat_0,g1_lon_1), \ float32 T_CDC_GDS1_EATM(g1_lat_0,g1_lon_1), float32 L_CDC_GDS1_LCY(g1_lat_0,g1_lon_1), \ float32 M_CDC_GDS1_MCY(g1_lat_0,g1_lon_1), float32 H_CDC_GDS1_HCY(g1_lat_0,g1_lon_1), \ float32 g1_lat_0(g1_lat_0), float32 g1_lon_1(g1_lon_1)'''print(fn.grib_source)'''FCST_d01_2018-07-27_12_00_00.grib'''
print(fn.variables.keys())'''odict_keys(['TMP_GDS1_HTGL', 'VIS_GDS1_SFC', 'R_H_GDS1_HTGL', 'A_PCP_GDS1_SFC_acc12h', 'NCPCP_GDS1_SFC_acc12h', 'ACPCP_GDS1_SFC_acc12h', 'T_CDC_GDS1_EATM', 'L_CDC_GDS1_LCY', 'M_CDC_GDS1_MCY', 'H_CDC_GDS1_HCY', 'g1_lat_0', 'g1_lon_1'])'''
tmp = fn.variables['TMP_GDS1_HTGL']print(tmp)"""float32 TMP_GDS1_HTGL(g1_lat_0, g1_lon_1) initial_time: 07/27/2018 (00:00) forecast_time_units: hours forecast_time: 12 level: 2 model: North Pacific Hurricane Wave Model parameter_number: 11 parameter_table_version: 2 gds_grid_type: 1 level_indicator: 105 _FillValue: 1e+20 units: K long_name: Temperature center: US National Weather Service - NCEP (WMC)unlimited dimensions: current shape = (921, 804)filling on"""
print(tmp.initial_time)'''07/27/2018 (00:00)'''
print(tmp[:]) print(tmp[:].data)
outfile = "test_out.nc"fw=nc.Dataset(outfile, 'w',format="NETCDF3_CLASSIC")
fw.MET_version = "V5.2"fw.model = "WRF"
nx = ny = 701fw.createDimension('lat',ny)fw.createDimension('lon',nx)
latitudes = fw.createVariable('lat', 'f4', ('lat',"lon",))longitudes = fw.createVariable('lon', 'f4', ("lat",'lon',))ftmax = fw.createVariable("TMAX","f4",("lat","lon",),fill_value=-9999)
latitudes.long_name="longitude"latitudes.units = "degrees_east"latitudes.standard_name="longitude"
longitudes.long_name="latitude"longitudes.units = "degrees_north"longitudes.standard_name="latitude"
ftmax.setncattr("name", "TMAX") ftmax.long_name = "Maximum temperature"ftmax.units = "K"ftmax.level = "Z2"ftmax.grib_code = 15
latitudes[:] = lats[:]longitudes[:] = lons[:]ftmax[:] = tmax[:]
fw.close()
import xarray as xr
ds = xr.open_dataset('ERA_HGT_P500_2019-08-02_00_00_00.nc')print(ds)"""Dimensions: (lat: 721, lon: 1440)Coordinates: * lat (lat) float32 90.0 89.75 89.5 89.25 ... -89.25 -89.5 -89.75 -90.0 * lon (lon) float32 0.0 0.25 0.5 0.75 1.0 ... 359.0 359.25 359.5 359.75Data variables: HGT (lat, lon) float32 ...Attributes: creation_date: Fri Aug 27 17:58:44 CST 2021 MET_version: V5.2 Nlon: 1440 grid_points Nlat: 721 grid_points delta_lon: 0.250000 degrees delta_lat: -0.250000 degrees lon_ll: 0.000000 degrees_east lat_ll: 90.000000 degrees_north Projection: LatLon"""
# 访问变量print(ds.HGT) # print(ds["HGT"])"""[1038240 values with dtype=float32]Coordinates: * lat (lat) float32 90.0 89.75 89.5 89.25 ... -89.25 -89.5 -89.75 -90.0 * lon (lon) float32 0.0 0.25 0.5 0.75 1.0 ... 359.0 359.25 359.5 359.75Attributes: valid_time_ut: 1564704000
valid_time: 20190802_000000 init_time_ut: 000000000.00000 init_time: 00000000_000000 grib_code: 7 units: gpm level: P500 long_name: Geopotential height name: HGT"""
ds.HGT.dims # dimensionsds.HGT.coords # Coordinate ds.HGT.attrs # attributes ds.HGT.data # data
ds.to_netcdf('output.nc')
往期推荐
