做ceemdan分解并写入nc
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 10 16:00:01 2017
@author: damtest
"""
import numpy as np
import pandas as pd
import netCDF4 as nc
from pyeemd import ceemdan
import os.path
fnc = nc.Dataset("/public/damtest/zy/data/station/cma.duan.sh.stn.73.year.mon.198001.201412.nc", "r")
sh=fnc.variables['sh']
stn=fnc.variables['stn']
ndims=sh.shape
nyear=ndims[0]
nmonth=ndims[1]
nstn=ndims[2]
data_in_memo=sh[:,:,:]
imfs=[]
for istn in range(nstn):
data=data_in_memo.reshape((nyear*nmonth,nstn))[:,istn]
imfs.append(ceemdan(data.filled(),ensemble_size=500))
sh_imfs=np.array(imfs)
try: ncfile.close() # just to be safe, make sure dataset is not already open.
except: pass
ncfile = nc.Dataset('tp_sh_stat_monthly_ceemdan_1979-2014.nc',mode='w',format='NETCDF4_CLASSIC')
stn_dim = ncfile.createDimension('stn', nstn) # latitude axis
imf_dim = ncfile.createDimension('imf', 8) # longitude axis
time_dim = ncfile.createDimension('time', None) # unlimited axis (can be appended to).
for dim in ncfile.dimensions.items():
print(dim)
ncfile.title="tp_sh_stat_monthly_ceemdan_1979-2014"
stn_out = ncfile.createVariable('stn', np.float32, ('stn',))
stn_out.units = 'n'
stn_out.long_name = 'stn number'
imf_out = ncfile.createVariable('imf', np.float32, ('imf',))
imf_out.units = 'n'
imf_out.long_name = 'imf '
time_out = ncfile.createVariable('time', np.float32, ('time',))
time_out.units = 'months since 1979-01-01'
time_out.long_name = 'time'
# Define a 3D variable to hold the data
sh_imfs_out = ncfile.createVariable('sh_imfs',np.float64,('stn','imf','time')) # note: unlimited dimension is leftmost
sh_imfs_out.units = 'w/m^2' # degrees Kelvin
sh_imfs_out.standard_name = 'sh over TP' # this is a CF standard name
print(sh_imfs_out)
print("-- Some pre-defined attributes for variable temp:")
print("temp.dimensions:", sh_imfs_out.dimensions)
print("temp.shape:", sh_imfs_out.shape)
print("temp.dtype:", sh_imfs_out.dtype)
print("temp.ndim:", sh_imfs_out.ndim)
stn_out[:] = stn[:]
imf_out[:] = range(8)
time_out[:] = range(nyear*nmonth)
sh_imfs_out[:,:,:] = sh_imfs
print("-- Wrote data, temp.shape is now ", sh_imfs_out.shape)
ncfile.close()
print('Dataset is closed!')
