"""Map 1D experimental profiles to netcdf grid generated by hypernotoad.pro
"""
from __future__ import (absolute_import, division,
print_function, unicode_literals)
__all__ = ["map_1d4grid",
"map_pfile2grid",
"map_nc2grid",
]
__date__ = '02072018'
__version__ = '0.2.1'
__author__ = 'J.G. Chen'
__email__ = 'cjgls@pku.edu.cn'
import os
import re
from shutil import copyfile
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as si
from scipy import ndimage, integrate
import boutpy.boututils.pfile
from boutpy.boututils.fileio import file_import, save2nc
from boutpy.boututils.functions import get_yesno
from boutpy.boutdata.boutgrid import boutgrid
[docs]def map_1d4grid(x, y, gridfile, index=None, ext='const', kind='cubic',
sigma=8, mode='nearest', wtitle='', Er=None, vfilter=None,
units=None,
smooth='1st', w=None, sfactor=None, verbose=False,
check_positive=False, **kwargs):
"""Map 1D-profile to grid
1. basic rule: constant in y-direction
2. map :math:`E_r`:
- E_r = er*Rxy*Bpxy/(Rxy*Bpxy)_omp,
- E_r = omgeb*Rxy*Bpxy
3. leg regions
- SOL region: same as the nearest line which has core part
- private flux region: constant values which are same as separatrix
.. warning::
The smooth method(`smooth=gaussian`: gaussian_filter1d) will slightly
make the gradient at boundary closer to 0.
Parameters
----------
x : (N,) array_like
1-D array of independent input data. Must be increasing.
y : (N,) array_like
1-D array of dependent input data, of the same length as x.
psin : (M,) array_like
1-D array of independent data.
index : int array_like or int, optional, default: None
ignoring data points at `index` when fitting data. if `index` is
a integer, then ignoring `index` data points at both 'x' boundary
for smoothing purposes.
smooth : ['gaussian', 'spline', '1st', '2nd'], optional, default: '1st'
smooth method if not `None`
- 'gaussian': scipy.ndimage.gaussian_filter1d
- 'spline': scipy.interpolate.UnivariantSpline
- '1st': smooth 1st order derivative
- '2nd': smooth 2nd order derivative
w : (M,) array_like, optional, default: None
weight for smoothing if 'spline' method used
sfactor : float, optional, default: None
smooth_factor for 'spline' method
ext : float, (float, float) or 'const', optional, default: 'const'
Controls the extrapolation mode for elements not in the interval
defined by the knot sequence.
- float : slope for the profiles outside the interval
- (float, float): slopes for both boundary
- 'const': constant value outside the interval
sigma : scalar, default: 8
standard deviation for Gaussian kernel
mode : ['reflect', 'constant', 'nearest', 'mirror', 'wrap'], optional
The `mode` parameter determines how the array borders are
handled, where `cval` is the value when mode is equal to
'constant'. Default is 'nearest'.
wtitle : str, optional, default: ''.
addtional title string for plot window
units : str, optional, default: None
additional info for the profile: units. It is shown in the fig title.
TODO: convert units in here?
Er : ['er', 'omgeb'], optional, default: None
- ``er``: map 'er' to 'E_r', E_r = er*Rxy*Bpxy/(Rxy*Bpxy)_omp
- ``omgeb``: map 'omgeb' to 'E_r', E_r = omgeb*Rxy*Bpxy
- ``None``: constant in y-direction.
vfilter : int|float, optional, default: None
Set the data where are less than `vfilter` to the one closest to
`vfilter`.
Not used by default.
It's recommended for electron temperature, especially when the Spitzer
resistivity, which is propotional to Te^{-1.5}, is used in the code.
check_positive : boel, optional, default: False
check that the profile only contains positive value. Raise ValueError
if it contains nagetive value or 0.
verbose : bool, optional, default: False
output more details
Other keywords arguments are passed to smooth method.
Returns
-------
data_2d : 2D-array_like
2D data after mapping 1D to grid
See Also
--------
scipy.ndimage.gaussian_filter1d,
scipy.ndimage.gaussian_filter,
scipy.interpolate.UnivariateSpline
"""
grid = boutgrid(gridfile)
# get psin at outer mid-plane
psin, yind_imp, yind_omp = grid.get_psin(
yind='omp', index=True, verbose=False)
ind_min = np.abs(x - psin[0]).argmin()
ind_max = np.abs(x - psin[-1]).argmin()
# convert unit to Latex format
if units:
# ^20 --> ^{20}
pat = re.compile(r'\^(-*\d+)')
for i in pat.findall(units):
units = units.replace('^{}'.format(i), '^{{{}}}'.format(i))
# _20 --> _{20}
pat = re.compile(r'_(-*\d+)')
for i in pat.findall(units):
units = units.replace('_{}'.format(i), '_{{{}}}'.format(i))
plt.figure('oned4grid-'+wtitle)
plt.clf()
fig, ax = plt.subplots(2, 2, sharex=True, num='oned4grid-'+wtitle)
ax = ax.flatten()
plt.suptitle(wtitle +
("(${}$)".format(units) if units is not None else ''))
ax[0].plot(x[ind_min:ind_max], y[ind_min:ind_max], 'k*-',
markevery=5, label='original')
grad_1st_tmp = np.gradient(y[ind_min:ind_max], x[ind_min:ind_max])
grad_2nd_tmp = np.gradient(grad_1st_tmp, x[ind_min:ind_max])
ax[1].plot(x[ind_min:ind_max], grad_1st_tmp, 'k*-',
markevery=5, label='original')
ax[2].plot(x[ind_min:ind_max], grad_2nd_tmp, 'k*-',
markevery=5, label='original')
# extend profile to outside the interval
data_1d = si.interp1d(x, y, kind='cubic', bounds_error=False,
fill_value=(y[0], y[-1]))(psin)
data_orignal = data_1d.copy()
# data_original[xind_l:xind_r+1] is experimental data
xind_l = 0
xind_r = len(psin) - 1
if x[0] > psin[0] or x[-1] < psin[-1]:
if x[0] > psin[0]:
# left boundary needs extension
xind_l = np.where(psin >= x[0])[0][0]
if x[-1] < psin[-1]:
# right boundary needs extension
xind_r = np.where(psin <= x[-1])[0][-1]
if isinstance(ext, basestring) and ext == 'const':
# constant outside the boundary
print("NOTE: set constant value outside the [{}, {}]!!!!!"
.format(x[0], x[-1]))
else:
if isinstance(ext, float):
slope_l = slope_r = ext
elif isinstance(ext, (list, tuple, np.ndarray)) and len(ext) == 2:
slope_l, slope_r = ext
else:
raise ValueError(
"option `ext`: a float value, 2-elements tuple of float "
"or 'const' required")
if xind_l != 0:
data_1d[:xind_l] = (slope_l * (psin[:xind_l] - psin[xind_l]) +
data_1d[xind_l])
if xind_r != len(psin) - 1:
data_1d[xind_r:] = (slope_r * (psin[xind_r:] - psin[xind_r]) +
data_1d[xind_r])
# using weight to remove some data
ind_keep = np.ones_like(psin)
xind_lpan = xind_rpan = None
if isinstance(index, (list, tuple, np.ndarray)):
ind_keep[np.asarray(index)] = 0
elif isinstance(index, int):
if index > psin.size / 8:
raise ValueError("option `index`: Too Large!")
if xind_l != 0:
xind_lpan = [max(xind_l-index, 1), min(xind_l+index, psin.size)]
ind_keep[xind_lpan[0]:xind_lpan[1]] = 0
if xind_r != len(psin) - 1:
xind_rpan = [max(xind_r-index, 1), min(xind_r+index, psin.size)]
ind_keep[xind_rpan[0]:xind_rpan[1]] = 0
if verbose:
print("[xind_l, xind_r] = [{}, {}]".format(xind_l, xind_r))
print("xind_lpan: ", xind_lpan)
print("xind_rpan: ", xind_rpan)
if (ind_keep == 0).any():
# some data points are removed
data_1d = si.interp1d(psin[np.where(ind_keep)],
data_1d[np.where(ind_keep)],
kind=kind)(psin)
if (isinstance(index, (tuple, list, np.ndarray))
and (np.diff(index) == 1).all()):
ax[0].axvspan(psin[index[0]], psin[index[-1]], alpha=0.5)
ax[1].axvspan(psin[index[0]], psin[index[-1]], alpha=0.5)
ax[2].axvspan(psin[index[0]], psin[index[-1]], alpha=0.5)
ax[3].axvspan(psin[index[0]], psin[index[-1]], alpha=0.5)
elif isinstance(index, int):
if xind_lpan is not None:
ax[0].axvspan(psin[xind_lpan[0]], psin[xind_lpan[-1]-1],
alpha=0.5)
ax[1].axvspan(psin[xind_lpan[0]], psin[xind_lpan[-1]-1],
alpha=0.5)
ax[2].axvspan(psin[xind_lpan[0]], psin[xind_lpan[-1]-1],
alpha=0.5)
ax[3].axvspan(psin[xind_lpan[0]], psin[xind_lpan[-1]-1],
alpha=0.5)
if xind_rpan is not None:
ax[0].axvspan(psin[xind_rpan[0]], psin[xind_rpan[-1]-1],
alpha=0.5)
ax[1].axvspan(psin[xind_rpan[0]], psin[xind_rpan[-1]-1],
alpha=0.5)
ax[2].axvspan(psin[xind_rpan[0]], psin[xind_rpan[-1]-1],
alpha=0.5)
ax[3].axvspan(psin[xind_rpan[0]], psin[xind_rpan[-1]-1],
alpha=0.5)
ax[0].plot(psin, data_1d, 'r+-', label='interp1d', markevery=7)
ax[1].plot(psin, np.gradient(data_1d, psin),
'r+-', label='interp1d', markevery=7)
ax[2].plot(psin, np.gradient(np.gradient(data_1d, psin), psin),
'r+-', label='interp1d', markevery=7)
ax[3].plot(psin[xind_l:xind_r+1],
((data_1d-data_orignal)/data_orignal)[xind_l:xind_r+1]*100,
'r+-', label='interp1d', markevery=7)
if isinstance(vfilter, (float, int)) and (data_1d < vfilter).any():
print("WARNING: set data where is less than {0} to {0}"
.format(vfilter))
ind = np.abs(data_1d - vfilter).argmin()
data_1d[np.where(data_1d < vfilter)] = data_1d[ind]
# weight for 'scipy.interpolate.UnivariateSpline()'
if w is None:
w = ind_keep
# smooth data, same psin range
grad_1st = np.gradient(data_1d, psin)
grad_2nd = np.gradient(grad_1st, psin)
if smooth == 'gaussian':
print("WARNING: gradient at boundary may be changed a lot!")
data_1d = ndimage.gaussian_filter1d(data_1d, sigma=sigma,
mode=mode, **kwargs)
elif smooth == 'spline':
spl = si.UnivariateSpline(psin, data_1d, w=w, **kwargs)
if sfactor:
spl.set_smoothing_factor(sfactor)
data_1d = spl(psin)
elif smooth == '1st':
# smooth 1st order derivative
s_grad_1st = ndimage.gaussian_filter1d(grad_1st, sigma=sigma,
mode=mode, **kwargs)
int1 = integrate.cumtrapz(s_grad_1st, psin, initial=0)
data_1d = int1 + (data_1d - int1).mean()
elif smooth == '2nd':
s_grad_2nd = ndimage.gaussian_filter1d(grad_2nd, sigma=sigma,
mode=mode, **kwargs)
int2_tmp = integrate.cumtrapz(s_grad_2nd, psin, initial=0)
int2_tmp += (grad_1st - int2_tmp).mean()
int2 = integrate.cumtrapz(int2_tmp, psin, initial=0)
data_1d = int2 + (data_1d - int2).mean()
if smooth:
ax[0].plot(psin, data_1d, 'b-', label=smooth)
ax[1].plot(psin, np.gradient(data_1d, psin), 'b-', label=smooth)
ax[2].plot(psin, np.gradient(np.gradient(data_1d, psin), psin),
'b-', label=smooth)
ax[3].plot(psin[xind_l:xind_r+1],
((data_1d-data_orignal)/data_orignal)[xind_l:xind_r+1]*100,
'b-', label=smooth)
ax[0].legend()
ax[-1].set_xlabel('$\psi_n$')
ax[-2].set_xlabel('$\psi_n$')
# ax[1].text(0.8, 0.8, 'df/d$\psi_n$', transform=ax[1].transAxes)
ax[1].legend([], title='df/d$\psi_n$')
ax[2].legend([], title='$d^2f/d\psi_n^2$')
ax[3].legend([], title='rel tol(%)')
if check_positive and (data_1d <= 0).any():
raise ValueError("negative value or 0 in profile!")
nx, ny = grid['nx'], grid['ny']
# psin[xsep] is larger than but nearest to 1.0??
xsep, xsep2 = grid.get('ixseps1', nx), grid.get('ixseps2', nx)
xsep = nx - 1 if xsep < 0 or xsep >= nx else xsep
xsep2 = nx - 1 if xsep2 < 0 or xsep2 >= nx else xsep2
assert xsep <= xsep2
nsol1 = nx - xsep
nsol2 = nx - xsep2
ny_inner = grid['ny_inner']
Rxy = grid['Rxy']
Bpxy = grid['Bpxy']
jyseps1_1 = grid['jyseps1_1']
jyseps2_1 = grid['jyseps2_1']
jyseps1_2 = grid['jyseps1_2']
jyseps2_2 = grid['jyseps2_2']
# initiate variable
data_2d = np.zeros([nx, ny])
# split yind
yind_leg1_1 = np.r_[0: jyseps1_1 + 1]
yind_leg2_2 = np.r_[jyseps2_2 + 1: ny]
if jyseps2_1 == jyseps1_2:
yind_leg2_1 = np.array([], dtype=int)
yind_leg1_2 = np.array([], dtype=int)
else:
yind_leg2_1 = np.r_[jyseps2_1 + 1: ny_inner]
yind_leg1_2 = np.r_[ny_inner: jyseps1_2 + 1]
yind_core = np.r_[jyseps1_1 + 1: jyseps2_1 + 1,
jyseps1_2 + 1: jyseps2_2 + 1]
yind_legd = np.r_[yind_leg1_1, yind_leg2_2]
yind_legu = np.r_[yind_leg2_1, yind_leg1_2]
yind_leg = np.r_[yind_legd, yind_legu]
if verbose >= 2:
print("yind for each region\n" + "=" * 18)
# print("yind_core: ", yind_core)
print("yind_legd: ", yind_legd)
print("yind_legu: ", yind_legu)
# print("yind_leg : ", yind_leg)
assert (len(yind_core) + len(yind_leg) == ny)
assert (len(yind_legd) + len(yind_legu) == len(yind_leg))
# start mapping
# -------------
# region 1: core + SOL
# constant in y-direction, except for E_r
if Er is None:
data_2d[:, yind_core] = data_1d.reshape(nx, 1)
elif Er == 'er':
print("NOTE: map er to E_r\n"
" E_r = er*Rxy*Bpxy/(Rxy*Bpxy)_omp\n")
data_2d[:, yind_core] = (data_1d.reshape(nx, 1)
* Rxy[:, yind_core]
* Bpxy[:, yind_core]
/ Rxy[:, yind_omp].reshape(nx, 1)
/ Bpxy[:, yind_omp].reshape(nx, 1)
)
elif Er == 'omgeb':
print("NOTE: map omgeb to E_r\n"
" E_r = omgeb*Rxy*Bpxy\n")
data_2d[:, yind_core] = (data_1d.reshape(nx, 1)
* Rxy[:, yind_core]
* Bpxy[:, yind_core])
else:
raise ValueError("ERROR: unrecognized value for option `Er`!!!")
# region 2: SOL regions of leg
# constant value at leg SOL region in y-direction
# set to the value nearest to region 1 respectively.
if Er is None:
data_2d[xsep:, yind_legd] = data_1d[xsep:].reshape(nsol1, 1)
data_2d[xsep2:, yind_legu] = data_1d[xsep2:].reshape(nsol2, 1)
else:
# E_r:
data_2d[xsep:, yind_leg1_1] = (
data_2d[xsep:, jyseps1_1 + 1].reshape(nsol1, 1))
data_2d[xsep:, yind_leg2_2] = (
data_2d[xsep:, jyseps2_2].reshape(nsol1, 1))
data_2d[xsep2:, yind_leg2_1] = (
data_2d[xsep2:, jyseps2_1].reshape(nsol2, 1))
data_2d[xsep2:, yind_leg1_2] = (
data_2d[xsep2:, jyseps1_2 + 1].reshape(nsol2, 1))
# region 3: private flux region
# constant: same value as separatrix
data_2d[:xsep, yind_legd] = data_2d[xsep, yind_legd]
data_2d[:xsep2, yind_legu] = data_2d[xsep2, yind_legu]
return data_2d
[docs]def map_pfile2grid(gridfile, pfile, output=None, Er='omgeb', verbose=False,
**kwargs):
"""Map pfile to netcdf grid generated by hypnotoad.pro
Parameters
----------
gridfile : str
path to gridfile
pfile : str
path to pfile
output : str, None, default: None
output data to file `output` which is a copy of `gridfile`,
if `output` is None, the suffix of filename is '_exp_constSOL_er'
(Er='er') or '_exp_constSOL_omgeb'(Er='omgeb'),
Er : ['er'|'omgeb'], default: 'omgeb'
- ``er``: map 'er' to 'E_r', E_r = er*Rxy*Bpxy/(Rxy*Bpxy)_omp
- ``omgeb``: map 'omgeb' to 'E_r', E_r = omgeb*Rxy*Bpxy
verbose : bool, optional, default: False
output more information.
Other keywords arguments are passed to map_1d4grid()
Returns
-------
None
"""
# load original data
grid = boutgrid(gridfile)
print("=" * 50)
psin, yind_imp, yind_omp = grid.get_psin(yind='omp', index=True)
xsep = np.where(psin >= 1.0)[0][0]
if grid.topology():
# with xpoint
assert xsep == grid['ixseps1']
else:
# no xpoint
assert xsep == psin.size
# interpolate all profiles in pfile to same psin
# try:
# pdata = file_import(pfile)
# except:
pfile_obj = boutpy.boututils.pfile(pfile)
pdata = pfile_obj.interpolate()
units = pfile_obj.units
profiles = {}
name = {'ti': 'Tiexp',
'te': 'Teexp',
'ni': 'Niexp',
'ne': 'Neexp',
'er': 'E_r',
'pressure_s': 'pressure_s'
}
print("smooth around separatrix area ...")
for i in ['ti', 'ni', 'ne']:
print("===> mapping {} ...".format(i))
profiles[i] = map_1d4grid(
pdata['psin'], pdata[i], grid,
wtitle=i, units=units[i], **kwargs)
# map te, filter out < 0.01 keV
if verbose >= 1:
print("===> mapping {} ...".format('te'))
profiles['te'] = map_1d4grid(
pdata['psin'], pdata['te'], grid, wtitle='te',
vfilter=0.05, units=units['te'], **kwargs)
if verbose >= 1:
print("===> mapping {} ...".format(Er))
profiles['er'] = map_1d4grid(
pdata['psin'], pdata[Er], grid,
wtitle=Er, Er=Er,
units=units[Er], **kwargs)
if verbose >= 1:
print("===> mapping {} ...".format('pressure'))
profiles['pressure_s'] = map_1d4grid(
psin, grid['pressure'][:, yind_omp],
grid, wtitle='pressure', units='Pa', **kwargs)
if verbose >= 1:
print("===> Unit conversion ...\n"
" er: kV/m --> V/m\n"
" ti, te: keV --> eV")
profiles['er'] *= 1000.0
profiles['ti'] *= 1000.0
profiles['te'] *= 1000.0
if isinstance(output, basestring):
output = output if output.endswith('.nc') else (output + '.nc')
else:
output = gridfile[:-3] + '_exp_constSOL_' + Er + '.nc'
if os.path.exists(output):
if not get_yesno('gridfile:\n {}\n'.format(output) +
'Already exists, Overwrite the file?'):
raise IOError("{} already exists!!!!".format(output))
# saving data to grid ...
print("save data to {} ...".format(output))
copyfile(gridfile, output)
save2nc(output, 'a', **{name[i]: profiles[i] for i in profiles.keys()})
[docs]def map_nc2grid(gridfile, ncfile, output=None, Er='omgeb', verbose=False,
**kwargs):
"""Map 1D profiles in `ncfile` to netcdf grid generated by hypnotoad.pro
Parameters
----------
gridfile : str
path to gridfile
ncfile : str
path to ncfile, which includes 1D profiles.
output : str, None, default: None
output data to file `output` which is a copy of `gridfile`,
if `output` is None, the suffix of filename is '_exp_constSOL_er'
(Er='er') or '_exp_constSOL_omgeb'(Er='omgeb'),
Er : ['er'|'omgeb'], default: 'omgeb'
'er': map 'er' to 'E_r', E_r = er*Rxy*Bpxy/(Rxy*Bpxy)_omp
'omgeb': map 'omgeb' to 'E_r', E_r = omgeb*Rxy*Bpxy
Other keywords arguments are passed to map_1d4grid()
Returns
-------
None
"""
# load original data
grid = boutgrid(gridfile)
# check topology of grid
nxpoint = grid.topology()
print("=" * 50)
psin, yind_imp, yind_omp = grid.get_psin(yind='omp', index=True)
xsep = np.where(psin >= 1.0)[0][0]
assert xsep == grid['ixseps1']
# interpolate all profiles in pfile to same psin
pdata = file_import(ncfile)
if (Er == 'omgeb') and (not pdata.has_key('omgeb')):
print("WARNING: `omgeb' not in {}".format(ncfile))
print(" Mapping Er using er")
Er = 'er'
if pdata['psin'][0] > psin[0]:
raise ValueError("psin range error")
units = {'ti': 'eV',
'te': 'eV',
'ne': '10^{20}/m^3',
'nee': '10^{20}/m^3',
'ni': '10^{20}/m^3',
'p0': 'Pa',
'er': 'V/m',
'omgeb': 'Rad/s',
}
profiles = {}
name = {'ti': 'Tiexp', # 'name in pdata': 'name in grid'
'te': 'Teexp',
'ni': 'Niexp',
'ne': 'Neexp',
'nee': 'Neexp',
'er': 'E_r',
'p0': 'pressure_s',
}
if verbose >= 1:
print("smooth around separatrix area ...")
for i in ['ti', 'ni', 'ne', 'nee', 'p0']:
if i not in pdata.keys():
continue
print("===> mapping {} ...".format(i))
profiles[i] = map_1d4grid(pdata['psin'], pdata[i], grid,
wtitle=i, units=units[i],
**kwargs)
# map te, filter out < 0.01
if verbose >= 1:
print("===> mapping {} ...".format('te'))
profiles['te'] = map_1d4grid(pdata['psin'], pdata['te'], grid,
wtitle='te', vfilter=0.01,
units=units['te'], **kwargs)
if verbose >= 1:
print("===> mapping {} ...".format(Er))
profiles['er'] = map_1d4grid(
pdata['psin'], pdata[Er], grid,
wtitle=Er, Er=Er,
units=units[Er], **kwargs)
if verbose >= 1:
print("===> mapping {} ...".format('pressure'))
if not profiles.has_key('p0'):
profiles['p0'] = map_1d4grid(
psin, grid['pressure'][:, yind_omp], grid, wtitle='pressure',
units='Pa', **kwargs)
if verbose >= 1:
print("===> Unit conversion ...")
if get_yesno("er: kV/m --> V/m?"):
profiles['er'] *= 1000.0
plt.figure('oned4grid-er')
# the plot use original data, so change to olddata
plt.title('er(kV/m)')
if get_yesno("omgeb: kRad/s --> Rad/s?"):
profiles['er'] *= 1000.0
plt.figure('oned4grid-omgeb')
plt.title('omgeb(kRad/s)')
if get_yesno("ti, te: keV --> eV?"):
profiles['ti'] *= 1000.0
profiles['te'] *= 1000.0
plt.figure('oned4grid-ti')
plt.title('ti(keV)')
if get_yesno("p0: kPa --> Pa?"):
profiles['p0'] *= 1000.0
plt.figure('oned4grid-pessure')
plt.title('pressure(kPa)')
if get_yesno("ne, ni: 10^19/m^3 --> 10^20/m^3?"):
profiles['nee'] /= 10.0
profiles['ni'] /= 10.0
plt.figure('oned4grid-nee')
plt.title('nee($10^{19}/m^3$)')
plt.figure('oned4grid-ni')
plt.title('ni($10^{19}/m^3$)')
if isinstance(output, basestring):
output = output if output.endswith('.nc') else (output + '.nc')
else:
if pdata['psin'][-1] < psin[-1]:
output = gridfile[:-3] + '_exp_constSOL_' + Er + '.nc'
else:
output = gridfile[:-3] + '_exp_' + Er + '.nc'
if os.path.exists(output):
if not get_yesno('gridfile:\n {}\n'.format(output) +
'Already exists, Overwrite the file?'):
raise IOError("{} already exists!!!".format(output))
# saving data to grid ...
print("save data to {} ...".format(output))
copyfile(gridfile, output)
save2nc(output, 'a', **{name[i]: profiles[i] for i in profiles.keys()})
