python 创建 读写 nc文件 demo

haysen

    做了两个demo，一个是用netcdf4模块，创建、读 nc文件，并保存为panel，一个是用xarray模块读和创建，基本都是各个模块说明书里的例子，自己做了demo自己备用，也分享给大家
https://github.com/haysengithub/python_creat_or_read_netcdf_demo.git

1. import time
   
2. import numpy as np
   
3. import pandas as pd
   
4. from netCDF4 import Dataset
   
5. 
   
6. class pyton_netcdf(object):
   
7.     def __init__(self):
   
8.         pass
   
9. 
   
10. 
    
11.     def creat_netcdf_by_netcdf4(self,file_name_nc):
    
12.         ntimes = 5
    
13.         nlevels = 10
    
14.         nlats = 73
    
15.         nlons = 144
    
16. 
    
17. 
    
18.         rootgrp = Dataset(file_name_nc, "w", format="NETCDF4")
    
19.         print(rootgrp.data_model)
    
20. 
    
21.         rootgrp.description = "bogus example script"
    
22.         rootgrp.history     = "Created"+time.ctime(time.time())
    
23.         rootgrp.source      = "netDF4 pytohn module tutorial"
    
24. 
    
25.         time_dim  = rootgrp.createDimension(dimname="time",  size=None)
    
26.         level_dim = rootgrp.createDimension(dimname="level", size=None)
    
27.         lat_dim   = rootgrp.createDimension(dimname="lat",   size=nlats)
    
28.         lon_dim   = rootgrp.createDimension(dimname="lon",   size=nlons)
    
29. 
    
30. 
    
31.         times = rootgrp.createVariable(varname="time", datatype="f8", dimensions=("time",))
    
32.         times.units = "hours since 0001-01-01 00:00:00.0"
    
33.         times.calendar = "gregorian"
    
34. 
    
35.         levels = rootgrp.createVariable(varname="level", datatype="i4", dimensions=("level",))
    
36.         levels.units = "hPa"
    
37. 
    
38.         lats = rootgrp.createVariable(varname="lat", datatype="f8", dimensions=("lat",))
    
39.         lats.units = "degrees north"
    
40. 
    
41.         lons = rootgrp.createVariable(varname="lon", datatype="f8", dimensions=("lon",))
    
42.         lons.units = "degrees east"
    
43. 
    
44.         temp = rootgrp.createVariable(varname="temp", datatype="f8", dimensions=("time", "level", "lat", "lon",))
    
45.         #压缩
    
46.         #temp = rootgrp.createVariable(varname="temp", datatype="f8", ("temp",),zlib=True, least_significant_digit=3)
    
47.         temp.units ="mm"
    
48. 
    
49. 
    
50. 
    
51.         from datetime import datetime, timedelta
    
52.         from netCDF4 import date2num
    
53.         times_data = [datetime(2001,3,1)+n*timedelta(hours=12) for n in range(ntimes)]
    
54.         times[:]   = date2num(times_data, units=times.units, calendar = times.calendar)
    
55. 
    
56.         levels[:]  = np.arange(start=0, stop = nlevels, step=1)
    
57. 
    
58.         lats[:]    = np.arange(start=-90, stop=91, step=2.5)
    
59.         lons[:]    = np.arange(start=-180,stop=180,step=2.5)
    
60. 
    
61.         from numpy .random import uniform
    
62.         temp[:] = uniform(size=(ntimes,nlevels,nlats,nlons))
    
63. 
    
64.         rootgrp.close()
    
65.         print(file_name_nc)
    
66. 
    
67.     def read_netcdf_by_netcdf4(self,file_name_nc):
    
68.         rootgrp = Dataset(file_name_nc)
    
69. 
    
70.         #print(rootgrp.dimensions)
    
71.         print("dimensions")
    
72.         for dimobj in rootgrp.dimensions.values():
    
73.             #print(dimobj)
    
74.             if dimobj.isunlimited():
    
75.                 print("\t\t", getattr(dimobj, "name"), "=UNLIMITED; //(", getattr(dimobj, "size"), "currently)")
    
76.                 #也可以
    
77.                 #print("\t\t", dimobj.name, "=UNLIMITED; //(", dimobj.size, "currently)")
    
78.             else:
    
79.                 print("\t\t", dimobj.name, "=", dimobj.size)
    
80. 
    
81.         #print(rootgrp.variables)
    
82.         print("variables:")
    
83.         for varobj in rootgrp.variables.values():
    
84.             print("\t\t", varobj.dtype, varobj.name, varobj.dimensions, ":")
    
85.             for var_cattrs_name in varobj.ncattrs():
    
86.                 print("\t\t"*2, varobj.name, ":", var_cattrs_name, "=", getattr(varobj,var_cattrs_name))
    
87. 
    
88.         print("//global attributes:")
    
89.         for global_attr in rootgrp.ncattrs():
    
90.             print("\t\t", global_attr, "=", getattr(rootgrp,global_attr))
    
91.         print("^"*20)
    
92. 
    
93.         # 直接读的例子
    
94.         print(rootgrp.description)
    
95.         print(rootgrp.dimensions["level"].name)
    
96.         print(rootgrp.dimensions["level"])
    
97.         print(rootgrp.variables["time"].units)
    
98.         print(rootgrp.variables["level"][:])
    
99.         print(rootgrp.variables["level"].shape)
    
100.         temp_nc_maskarray = rootgrp.variables["temp"][:]
     
101.         print(temp_nc_maskarray)
     
102.         print(temp_nc_maskarray.data)
     
103.         print(temp_nc_maskarray.mask)
     
104.         print(temp_nc_maskarray.fill_value)
     
105. 
     
106. 
     
107.     def netcdf_to_panel_by_netcdf4(self,file_name_nc):
     
108.         rootgrp = Dataset(file_name_nc)
     
109. 
     
110. 
     
111.         temp_panel = pd.Panel(data=rootgrp.variables["temp"][:, 0, :, :].data,
     
112.                               items= rootgrp.variables["time"][:].data,
     
113.                               major_axis=rootgrp.variables["lat"][:].data,
     
114.                               minor_axis=rootgrp.variables["lon"][:].data)
     
115.         print(temp_panel)
     
116. 
     
117.         from netCDF4 import num2date
     
118. 
     
119.         times = rootgrp.variables["time"]
     
120. 
     
121.         times_data = num2date(times[:], units=times.units, calendar=times.calendar)
     
122.         print(times_data)
     
123.         '''
     
124.         temp_panel4d = pd.Panel4D(data=rootgrp.variables["temp"][:, :, :].data,
     
125.                                 labels=pd.period_range(start=times_data[0].strftime("%Y%m%d%H%M"), end=None, periods=times.shape[0],freq="12h",name="time"),
     
126.                                 items=rootgrp.variables["level"][:].data,
     
127.                                 major_axis=rootgrp.variables["lat"][:].data,
     
128.                                 minor_axis=rootgrp.variables["lon"][:].data)
     
129.         print(temp_panel4d)
     
130.         '''
     
131. 
     
132.     def creat_netcdf_by_xarray(self, file_name_nc):
     
133.         #Creat an example dateset
     
134.         import numpy as np
     
135.         import pandas as pd
     
136.         import xarray as xr
     
137. 
     
138.         ntimes = 5
     
139.         nlevels = 10
     
140.         nlats = 73
     
141.         nlons = 144
     
142. 
     
143.         levels  = np.arange(start=0, stop = nlevels, step=1)
     
144. 
     
145.         lats    = np.arange(start=-90, stop=91, step=2.5)
     
146.         lons    = np.arange(start=-180,stop=180,step=2.5)
     
147. 
     
148.         from numpy .random import uniform
     
149.         temps = uniform(size=(ntimes,nlevels,nlats,nlons))
     
150. 
     
151.         times = pd.date_range("2001-03-01 00:00", periods = temps.shape[0],freq="12h")
     
152. 
     
153.         dims = ["time", "level", "lat", "lon"]
     
154. 
     
155.         temp_attr = dict(standard_name="air_potential_tempeature", uints="f")
     
156.         prcp_attr = dict(standard_name="convective_precipitation_flux", uints="mm")
     
157. 
     
158.         level_attr = dict(standard_name="level", uints="hPa")
     
159.         lat_attr   = dict(standard_name="lat", uints="degrees north")
     
160.         lon_attr   = dict(standard_name="lon", uints="degrees east")
     
161. 
     
162.         ds = xr.Dataset({
     
163.             "temperature":(dims,temps,temp_attr),
     
164.             "precipitation":(dims,temps,prcp_attr)},
     
165.             coords={
     
166.                 "time":times,
     
167.                 "level":(["level"],levels,level_attr),
     
168.                 "lat":(["lat"],lats,lat_attr),
     
169.                 "lon":(["lon"],lons,lon_attr)
     
170.                 #"reference_time":pd.Timestamp("2001-03-01 00:00")
     
171.                  }
     
172.             )
     
173. 
     
174.         ds.attrs=dict(description="bogus example script",
     
175.                     history="Created"+time.ctime(time.time()),
     
176.                     source="xarray python module tuorial",
     
177.                     author="some one")
     
178. 
     
179.         ds.to_netcdf(path=file_name_nc, mode="w", format="NETCDF4")
     
180. 
     
181.     def read_netcdf_to_xarray(self, file_name_nc):
     
182.         import xarray
     
183.         ds = xarray.open_dataset(file_name_nc)
     
184. 
     
185.         print("-"*20)
     
186. 
     
187.         #查看attrs，orderedDict，可以按照字典进行遍历
     
188.         print(ds.attrs)
     
189.         print("-"*20)
     
190. 
     
191.         #查看dims，SortedKeysDict，可以按照字典进行遍历
     
192.         print(ds.dims)
     
193.         print("-"*20)
     
194. 
     
195.         #查看coordinates 可以按照字典遍历，所谓coordinates，这里专指time level lat lon 这样的维度的数据
     
196.         print(ds.coords)
     
197.         print("-"*20)
     
198. 
     
199.         #coords_value 是一个DataArray
     
200.         for coords_value in ds.coords.values():
     
201.             print(coords_value)
     
202.             print("-"*5)
     
203.             print(coords_value.dims)
     
204.             print("-"*5)
     
205.             print(coords_value.values)
     
206.             print("-"*5)
     
207.             print(coords_value.coords)
     
208.             print("-"*5)
     
209.             print(coords_value.name)
     
210.             print("-"*5)
     
211.             print(coords_value.attrs)
     
212.             print("-"*5)
     
213.             print(coords_value.data)
     
214.         print("-"*20)
     
215. 
     
216.         #举个例子
     
217.         print(ds.coords["lon"].attrs)
     
218.         print(ds.coords["level"].values)
     
219. 
     
220.         print("-"*20)
     
221. 
     
222.         #查看数据，这个是除了维度信息后真正的数据，比如本例子中的温度和降雨量，也是个字典
     
223.         for data_value in ds.data_vars.values():
     
224.             print(data_value)
     
225.             print("-"*5)
     
226.             print(data_value.dims)
     
227.             print("-"*5)
     
228.             print(data_value.values)
     
229.             print("-"*5)
     
230.             print(data_value.coords)
     
231.             print("-"*5)
     
232.             print(data_value.name)
     
233.             print("-"*5)
     
234.             print(data_value.attrs)
     
235.             print("-"*5)
     
236.             print(data_value.data)
     
237. 
     
238.             #再举个例子
     
239.             print(ds.data_vars["precipitation"].values)  #同data
     
240.             print(ds.data_vars["precipitation"].attrs)   #参数
     
241.             print(ds.data_vars["precipitation"].name)    #名字
     
242.             print(ds.data_vars["precipitation"].data)    #真正的数据
     
243. 
     
244.             print(ds.data_vars["temperature"].values)
     
245.             print(ds.data_vars["temperature"].attrs)
     
246.             print(ds.data_vars["temperature"].name)
     
247.             print(ds.data_vars["temperature"].data)
     
248. 
     
249. 
     
250.         #这条命令类似于于ncdump -h
     
251.         print(ds.info())
     
252. 
     
253. 
     
254. 
     
255. 
     
256. 
     
257. 
     
258. 
     
259. 
     
260. 
     
261. 
     
262. 
     
263. if __name__ == "__main__":
     
264. 
     
265.     py_nc = pyton_netcdf()
     
266.     py_nc.creat_netcdf_by_netcdf4("netcdf4.nc")
     
267.     py_nc.read_netcdf_by_netcdf4("netcdf4.nc")
     
268.     py_nc.netcdf_to_panel_by_netcdf4("netcdf4.nc")
     
269. 
     
270.     py_nc.creat_netcdf_by_xarray("xarray.nc")
     
271.     py_nc.read_netcdf_to_xarray("xarray.nc")

复制代码

又是那隻貓

写的太复杂了，用函数会更好点

931404656

猫神可以给点例子+建议

haysen

又是那隻貓 发表于 2018-8-30 22:08
写的太复杂了，用函数会更好点

嗯，这个主要是给自己看的，写这样是很多用法是从help文件里拔出来的，怕自己忘了，干脆都写上，对我自己是一目了然，对别人可能会感觉冗余，稍后有空再出个精简版的，后面还会补充时间序列处理，比如小时转日数量，日转月啥的，不过我看NCL这些函数都有了，打算还是直接用ncl了

xusuperver

牛逼，问一下，哪里下载necp资料

尉鹏

python 可以用pynio

Benjamin_19

有点复杂，没看懂

正能量

良心帖子，非常好！！！

haysen

有位 冬青老先生 问我问题，我发现我没有发消息的权限（这个设定是脑子抽了么，别人可以给我发消息，我没有权限回人家）
就在这里回复了
你好，没有前后语境我也看不出来哦，目前看到的units 应该不对，这个要符合nc文件时间格式，比如我的例子里 "hours since 0001-01-01 00:00:00.0"   ，或者"days since 2017-01-01" 这样， 你只有一个hours是不对的吧。你的time数组我看不到。如果你用netcdf4这个包的话，你可以help(num2date) 查看它的用法

joannaxuxu

解决了我不会读nc文件全局属性的问题！{:eb502:}