python 写nc的示例
# coding=utf8
'''
creat by yushuai
time 2019 01 18
'''
import pdb
import sys
import os
import time
from netCDF4 import Dataset
import scipy.ndimage
import numpy as np
from datetime import datetime
from datetime import timedelta
def get_files(inpath):
files=os.listdir(inpath)
files.sort()
return files
def load_data(infile,element):
try:
nc = Dataset(infile, "r")
‘’‘
在读nc时候有时候 有时需要屏蔽mask和 放大缩小倍数
or
nc.set_auto_mask(False)
’‘’
if element=='total_cloud_cover':
# lat = nc.variables['LAT'][:]
# lon = nc.variables['LON'][:]
data = nc.variables[element][:]
else:
data=nc.variables[element][:][33, :]
nc.close()
except:
raise ValueError('Cant read this data')
return data
def chazhi_zoom(data,res):
'''
双线性插值将是order = 1
最临近插值的是order = 0
立方体是默认值（顺序= 3）
'''
return scipy.ndimage.zoom(data, res, order=1)
def creat_nc(filename,data,element_name,outpath):
lat = np.linspace(0,60.0,2401)
lon = np.linspace(70.0,140.0,2801)
# lon,lat=np.meshgrid(lon,lat)
if not os.path.exists(outpath+os.sep+filename):
da=Dataset(outpath+os.sep+filename,"w",format="NETCDF4")
da.createDimension("lat",lat.shape[0])
da.createDimension("lon",lon.shape[0])
t3=da.createVariable(element_name,"f8",("lat","lon"))
lonn=da.createVariable("lon","f8",("lon",))
latt=da.createVariable("lat","f8",("lat",))
da.variables[element_name][:]=data
da.variables["lat"][:]=lat
da.variables["lon"][:]=lon
#子属性
# t3.name="t3"
t3.level = element_name
t3.init_time = "20050807_000000"
t3.init_time_ut = "1123372800"
t3.valid_time = "20050807_120000"
t3.valid_time_ut = "1123416000"
t3.accum_time = "120000"
t3.setncattr('accum_time_sec', np.int32(43200) )
# t3.accum_time_sec = 1
latt.unit='degrees_north'
latt.long_name = "latitude"
latt.standard_name = "latitude"
lonn.unit='degrees_east'
lonn.long_name = "longitude"
lonn.standard_name = "longitude"
#全局属性
da.Projection = "LatLon" #"Lambert Conformal"
da.MET_version = "V8.0"
# da.MET_tool = "pcp_combine"
da.lat_ll = "0.0000 degrees_north"
da.lon_ll = "70.0000 degrees_east"
da.delta_lat = "0.0250 degrees"
da.delta_lon = "0.0250 degrees"
da.Nlat = "2401 grid_points"
da.Nlon = "2801 grid_points"
da.close()
def main1():
inpath=r'/g4/zhuzhi/kts/laps_data_3km'
outpath=r'/g4/zhuzhi/kts/yushuai/data/met_laps_3km'
date_s='20180202'
date_e='20180202'
date_s = datetime.strptime(date_s, '%Y%m%d')
date_e = datetime.strptime(date_e, '%Y%m%d')
while date_s <=date_e:
time= date_s.strftime('%Y%m%d')
opath=outpath+os.sep+time
if not os.path.exists(opath):
os.makedirs(opath)
path=inpath+os.sep+time+os.sep+"data"
files=get_files(path)
for file in files:
element=''
if file.startswith( 'u3' ) or file.startswith( 'v3' ) or file.startswith( 'lcp' ):
element=file.split("_")[0]
elif file.startswith( 'lcv' ):
element='total_cloud_cover'
elif file.startswith( 'lt1' ):
element=file.split("_")[1]
if element:
infile=path+os.sep+file
print infile,element
data=load_data(infile,element)
creat_nc(file,data,element,opath)
aday = timedelta(days=1)
date_s = date_s + aday
def main2():
inpath=r'/g4/zhuzhi/kts/laps_data_3km/5km'
outpath=r'/g4/zhuzhi/kts/yushuai/data/met_laps_3km'
files=get_files(inpath)
for file in files:
try:
element=''
if file.startswith('lh3') :
element=file.split("_")[0]
time=file.split("_")[1].split(".")[0][:-2]
# elif file.startswith('3d_rai'):
# element=file.split("_")[0]+'_'+file.split("_")[1]
# res=5
# time=file.split("_")[2].split(".")[0][:-2]
if element:
infile=inpath+os.sep+file
opath=outpath+os.sep+time
if not os.path.exists(opath):
os.makedirs(opath)
print infile,element
if not os.path.exists(opath+os.sep+file):
data=load_data(infile,element)
creat_nc(file,data,element,opath)
else:
pass
except:
pass
if __name__ == '__main__':
main1()
# print files
# data = load_u3_data()
# pdb.set_trace()
# creat_nc('u3_2018020100.nc',data,'u3')
# infile=r'E:\TEST\data\laps\lt1_t3_2018020104.nc'
# element='t3'
# res=2.5
# data=load_data(infile,element,res)
# file='lt1_t3_2018020104.nc'
# opath=r'E:\TEST\data\laps_new'
# creat_nc_new(file,data,element,opath)
