本文整理汇总了Python中xarray.open_dataarray函数的典型用法代码示例。如果您正苦于以下问题:Python open_dataarray函数的具体用法?Python open_dataarray怎么用?Python open_dataarray使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了open_dataarray函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的Python代码示例。
示例1: timeline_trend_count_SA
def timeline_trend_count_SA():
msg_folder = cnst.GRIDSAT
fname = 'aggs/gridsat_WA_-65_monthly_count_-40base_1000km2.nc'
fname2 = 'aggs/gridsat_WA_-40_monthly_count_-40base_1000km2.nc'
da = xr.open_dataarray(msg_folder + fname)
da2 = xr.open_dataarray(msg_folder + fname2)
#[25,33,-28,-10] , West[15,25,-26,-18]
da = da.sel(lat=slice(-25,-18), lon=slice(18, 22))# (lat=slice(-28,-10), lon=slice(25, 33))
da2 = da2.sel(lat=slice(-25,-18), lon=slice(18, 22)) #[25,33,-28,-10]
#da=da.sel(lat=slice(5,10))
#da[da==0]=np.nan
mean = da.mean(dim=['lat', 'lon'])
mean2 = da2.mean(dim=['lat', 'lon'])
#mean = mean[(mean['time.month']==8)]
f= plt.figure(figsize=(10,6))
for i in [12,1]:
bla = mean[(mean['time.month'] == i)]
bla.plot(label=str(i), marker='o')
plt.title('Average number of pixels <= -70C, SouthA 10-28S, 25-35E')
f = plt.figure(figsize=(10, 6))
for i in [12,1]:
bla2 = mean2[(mean2['time.month'] == i)]
bla2.plot(label=str(i), marker='o')
plt.title('Average number of pixels <= -40C, SouthA 10-28S, 25-35E')
plt.legend()
开发者ID:cornkle,项目名称:proj_CEH,代码行数:29,代码来源:gridsat_baseplots.py
示例2: run
def run(params):
start_time = datetime.now()
bin_width, filter_bandwidth, theta, shift, signal_field = params
# Get file paths
signal_dir = '/scratch/pkittiwi/fg1p/signal_map/bin{:.2f}/' \
'fbw{:.2f}/theta{:.1f}/shift{:d}' \
.format(bin_width, filter_bandwidth, theta, shift)
output_dir = '/scratch/pkittiwi/fg1p/stats_semi/signal/bin{:.2f}/' \
'fbw{:.2f}/theta{:.1f}/shift{:d}' \
.format(bin_width, filter_bandwidth, theta, shift)
signal_file = '{:s}/signal_map_bin{:.2f}_fbw{:.2f}_' \
'theta{:.1f}_shift{:d}_{:03d}.nc'\
.format(signal_dir, bin_width, filter_bandwidth,
theta, shift, signal_field)
output_file = '{:s}/stats_semi_signal_bin{:.2f}_fbw{:.2f}_' \
'theta{:.1f}_shift{:d}_{:03d}.nc' \
.format(output_dir, bin_width, filter_bandwidth,
theta, shift, signal_field)
mask_file = '/scratch/pkittiwi/fg1p/hera331_fov_mask.nc'
# Load data to memory and align coordinates
with xr.open_dataarray(signal_file) as da:
signal = da.load()
with xr.open_dataarray(mask_file) as da:
mask = da.load()
# Load one noise file to get coordinates.
noise = xr.open_dataarray(
'/scratch/pkittiwi/fg1p/noise_map/bin0.08/fbw8.00/theta90.0/shift0/'
'noise_map_bin0.08_fbw8.00_theta90.0_shift0_333.nc'
)
for key, values in noise.coords.items():
signal.coords[key] = values
mask.coords[key] = values
signal, noise, mask = xr.align(signal, noise, mask)
# Mask observation
signal = signal.where(mask == 1)
# Calculate statistic
out = get_stats(signal)
out.attrs = {'bin_width': bin_width, 'filter_bandwidth': filter_bandwidth,
'theta': theta, 'shift': shift}
os.makedirs(output_dir, exist_ok=True)
out.to_netcdf(output_file)
out.close()
print('Finish. signal_file = {:s}. output_file = {:s}. '
'Time spent {:.5f} sec.'
.format(signal_file, output_file,
(datetime.now() - start_time).total_seconds()))
开发者ID:piyanatk,项目名称:sim,代码行数:54,代码来源:calculate_stats_semi_mp.py
示例3: __init__
def __init__(self, instrument, ref_pb=None):
"""
Parameters
----------
:param instrument : Instrument
Instrument configuration
:param ref_pb : str, optional
name of the reference passband
"""
super().__init__(instrument, ref_pb)
I = self.instrument
self._spectra = sp = pd.read_hdf(resource_filename(__name__, join("data", "spectra.h5")), 'Z0')
self._tr_table = trt = xa.open_dataarray(resource_filename(__name__, join("data", "transmission.nc")))
self._tr_mean = trm = trt.mean(['airmass', 'pwv'])
self.extinction = interp1d(trm.wavelength, trm, bounds_error=False, fill_value=0.0)
self.wl = wl = sp.index.values
self.lte = lte = sp.columns.values
self._apply_extinction = True
# Dataframe indices
# -----------------
self.ipb = pd.Index(self.instrument.pb_names, name='passband')
self.iteff = pd.Index(lte, name='teff')
# Per-passband fluxes
# -------------------
self._compute_relative_flux_tables(0, ref_pb)
开发者ID:hpparvi,项目名称:PyTransit,代码行数:29,代码来源:contamination.py
示例4: composite
def composite(h):
pool = multiprocessing.Pool(processes=8)
file = '/users/global/cornkle/MCSfiles/blob_map_allscales_-50_JJAS_points_dominant.nc'
msg = xr.open_dataarray(file)
msg = msg[(msg['time.hour'] == 17) & (msg['time.minute'] == 0) & (
msg['time.year'] >= 2006) & (msg['time.year'] <= 2009) & (msg['time.month'] >= 6) ]
msg = msg.sel(lat=slice(10.5,17.5), lon=slice(-9.5,9.5))
res = pool.map(file_loop, msg)
pool.close()
# res = []
# for m in msg[0:50]:
# r = file_loop(m)
# res.append(r)
res = [x for x in res if x is not None]
scales = res
scales = [item for sublist in scales for item in sublist] # flatten list of lists
scales = np.concatenate(scales)
return scales
开发者ID:cornkle,项目名称:proj_CEH,代码行数:32,代码来源:surfaceScales_powerValues.py
示例5: saveDailyBlobs
def saveDailyBlobs():
"""
Converts hourly centre-point convective-core files to daily netcdf files so they can be saved with LSTA daily data
:return:
"""
msgfile = '/users/global/cornkle/MCSfiles/blob_map_allscales_-50_JJAS_points_dominant.nc'
msg = xr.open_dataarray(msgfile)
# def first_nozero(array_like, axis):
# array_like[array_like<16]= array_like[array_like<16]+24
# return np.nanmin(array_like,axis=axis)
msg.values[msg.values > 75] = np.nan
msg.values[msg.values == 0] = np.nan
for m in msg:
if m['time.hour'].values >= 16:
m.values[m > 0] = m['time.hour'].values
else:
m.values[m > 0] = m['time.hour'].values+24
### this is useful, it removes all pixels which got rain twice on a day
md = msg.resample('24H', base=16, dim='time', skipna=True, how='min')
md = md[(md['time.month'] >=6) & (md['time.month'] <=9)]
md.values[md.values>23] = md.values[md.values>23]-24
md.to_netcdf('/users/global/cornkle/MCSfiles/blob_map_allscales_-50_JJAS_points_dominant_daily.nc')
开发者ID:cornkle,项目名称:proj_CEH,代码行数:30,代码来源:saveLSTA.py
示例6: setUp
def setUp(self):
file = os.path.join(BASE_PATH, 'model', 'GFS_Global_0p25deg_20161219_0600.nc')
ds = SpatialDataset(NetCDFHandler(file),)
self.array = xr.open_dataarray(file)
self.grid = ds.get_grid(
'Maximum_temperature_height_above_ground_Mixed_intervals_Maximum',
data_array=self.array)
开发者ID:maestrotf,项目名称:pymepps,代码行数:7,代码来源:test_spatial.py
示例7: test_save_saves_also_grid
def test_save_saves_also_grid(self):
self.array.pp.grid = self.grid
self.array.pp.save('test.nc')
opened_array = xr.open_dataarray('test.nc')
grid_attrs = {attr[7:]: opened_array.attrs[attr]
for attr in opened_array.attrs if attr[:7] == 'ppgrid_'}
opened_grid = GridBuilder(grid_attrs).build_grid()
self.assertEqual(self.grid, opened_grid)
开发者ID:maestrotf,项目名称:pymepps,代码行数:8,代码来源:test_spatial.py
示例8: run_rebuild_iter_1mhz
def run_rebuild_iter_1mhz(args):
process = multiprocessing.current_process().pid
inf, ouf = args
print('pid: {:d} ; input file: {:s} ; output file: {:s}'
.format(process, inf, ouf))
w_in = xr.open_dataarray(inf)
w_out = rebuild_iter_1mhz(w_in)
w_out.to_netcdf(ouf)
开发者ID:piyanatk,项目名称:sim,代码行数:8,代码来源:resameple_filter.py
示例9: load
def load(cls, name, **kwargs):
Bx = xr.open_dataarray(name, **kwargs)
spl = cls(knots=Bx.knots, order=Bx.order,
bc=Bx.bc, dim=Bx.dim)
spl._coef = Bx
return spl
开发者ID:nbren12,项目名称:gnl,代码行数:9,代码来源:xspline.py
示例10: test_load_removes_grid_attrs
def test_load_removes_grid_attrs(self):
self.array.pp.grid = self.grid
self.array.pp.save('test.nc')
non_gridded_array = xr.open_dataarray('test.nc')
non_gridded_attrs = [attr for attr in non_gridded_array.attrs
if attr[:7] == 'ppgrid_']
self.assertTrue(non_gridded_attrs)
gridded_array = xr.DataArray.pp.load('test.nc')
gridded_attrs = [attr for attr in gridded_array.attrs
if attr[:7] == 'ppgrid_']
self.assertFalse(gridded_attrs)
开发者ID:maestrotf,项目名称:pymepps,代码行数:11,代码来源:test_spatial.py
示例11: regrid_simpler
def regrid_simpler(cmorph):
dummy = xr.open_dataset(constants.LSTA_TESTFILE)
cm = xr.open_dataarray(cmorph)
out = cmorph.replace('WA_', 'WA_onLSTA_')
cm_on_lst = dummy.salem.transform(cm)
enc = {'pr': {'complevel': 5, 'zlib': True}}
cm_on_lst.to_netcdf(out, encoding=enc, format='NETCDF4')
开发者ID:cornkle,项目名称:proj_CEH,代码行数:11,代码来源:saveCMORPH.py
示例12: run
def run(binnum):
process = current_process().pid
print('... P{:d}: applying filter {:s}'
.format(process, filter_files[binnum].split('/')[-1]))
filter_da = xr.open_dataarray(filter_files[binnum])
filter_array = filter_da.values
data_channels = filter_da.attrs['frequency_channels']
filter_bandwidth = filter_da.attrs['filter_bandwidth']
# Figure out FFT and filter normalization
# FFT normalization factor
x = filter_da.attrs['x']
y = filter_da.attrs['y']
f = filter_da.attrs['f']
dx = x[1] - x[0]
dy = y[1] - y[0]
df = f[1] - f[0]
u = filter_da.attrs['u']
v = filter_da.attrs['v']
e = filter_da.attrs['e']
du = u[1] - u[0]
dv = v[1] - v[0]
de = e[1] - e[0]
fft_norm = dx * dy * df
ifft_norm = du * dv * de * filter_array.size
# Filter normalization factor
filter_volume = np.sum(filter_array.size * du * dv * de)
filter_integral = np.sum(np.abs(filter_array) ** 2 * du * dv * de)
filter_norm = np.sqrt(filter_volume / filter_integral)
# Apply filter
filtered_data = apply_filter(
data_array[data_channels], filter_array,
fft_multiplier=fft_norm, ifft_multiplier=ifft_norm,
output_multiplier=filter_norm, apply_window_func=args.apply_window_func,
invert_filter=False
).real
out_da_attrs = filter_da.attrs
out_da_attrs.pop('x')
out_da_attrs.pop('y')
out_da_attrs.pop('f')
out_da_attrs['kx'] = filter_da.kx.values
out_da_attrs['ky'] = filter_da.ky.values
out_da_attrs['kz'] = filter_da.kz.values
out_da = xr.DataArray(
filtered_data, dims=['f', 'y', 'x'], coords={'f': f, 'y': y, 'x': x},
attrs=out_da_attrs
)
outfile = '{:s}/signal_cube_filtered_fbw{:.2f}MHz_{:03d}_bin{:03d}.nc'\
.format(args.output_directory, filter_bandwidth / 1e6,
field_num, binnum)
out_da.to_netcdf(outfile)
开发者ID:piyanatk,项目名称:sim,代码行数:53,代码来源:apply_filters_nostitch.py
示例13: composite
def composite():
pool = multiprocessing.Pool(processes=4)
file = constants.MCS_POINTS_DOM
msg = xr.open_dataarray(file)
msg = msg[ (msg['time.minute'] == 0) & (
msg['time.year'] >= 2006) & (msg['time.year'] <= 2010) & (msg['time.month'] >= 6)] #(msg['time.hour'] >= 17) &
msg = msg.sel(lat=slice(10.2, 17), lon=slice(-9.5, 9.5))
res = pool.map(file_loop, msg)
pool.close()
# for m in msg[0:10]:
# file_loop(m)
#
# return
res = [x for x in res if x is not None]
cell = []
surface = []
hour = []
for r in res:
cell.append(r[0])
surface.append(r[1])
hour.append(r[2])
pdb.set_trace()
cell = [item for sublist in cell for item in sublist] # flatten list of lists
surface = [item for sublist in surface for item in sublist] # flatten list of lists
hour = [item for sublist in hour for item in sublist] # flatten list of lists
cell = np.array(cell, dtype=float)
cell = cell[np.isfinite(surface)]
surface = np.array(surface, dtype=float)
surface = surface[np.isfinite(surface)]
hour = np.array(hour, dtype=float)
hour = hour[np.isfinite(surface)]
dic = {'cell': cell,
'surface': surface,
'hour' : hour }
pkl.dump(dic,
open("/users/global/cornkle/figs/LSTA-bullshit/scales/new/dominant_scales_save/scatter_scales.p", "wb"))
print('Successfully written scatter_scales save file')
开发者ID:cornkle,项目名称:proj_CEH,代码行数:50,代码来源:scatter_scales.py
示例14: blobs
def blobs():
#file = '/users/global/cornkle/MCSfiles/blob_map_30km_-67_JJAS_points.nc'
file = '/users/global/cornkle/MCSfiles/blob_map_allscales_-50_JJAS_points_dominant.nc'
fpath = '/users/global/cornkle/data/pythonWorkspace/proj_CEH/topo/gtopo_1min_afr.nc'
msg = xr.open_dataarray(file)
msg = msg.sel(lat=slice(10, 20), lon=slice(-10, 10))
msg = msg[ (msg['time.month'] >= 6 ) ]
msg = msg.where(msg > 6)
msg.values[msg.values>6] = 1
msg = msg.sum(dim='time')
map = msg.salem.get_map(cmap='viridis')
top = xr.open_dataarray(fpath)
f = plt.figure()
z = map.set_topography(top, relief_factor=1.4)
map.set_contour(z, levels=(200,400,600,800), cmap='Reds' )
map.set_data(msg)
map.visualize(title='Blobs and topo')
msg = msg.sum(dim='lon')
f = plt.figure()
msg.plot()
开发者ID:cornkle,项目名称:proj_CEH,代码行数:24,代码来源:test_lsta.py
示例15: get_previous_hours_msg
def get_previous_hours_msg(date, ehour, refhour):
# tdic = {18 : ('36 hours', '15 hours'),
# 19 : ('37 hours', '16 hours'),
# 20: ('38 hours', '17 hours'),
# 21: ('39 hours', '18 hours'),
# 22: ('40 hours', '19 hours'),
# 23: ('41 hours', '20 hours'),
# 0: ('42 hours', '21 hours'),
# 3: ('45 hours', '24 hours'),
# 6: ('48 hours', '27 hours')}
# before = pd.Timedelta(tdic[date.hour][0])
# before2 = pd.Timedelta(tdic[date.hour][1])
date = date.replace(hour=refhour)
if ehour > 0:
edate = date + pd.Timedelta(str(ehour) + ' hours')
else:
edate = date - pd.Timedelta(str(np.abs(ehour)) + ' hours')
#edate = edate.replace(hour=ehour)
t1 = edate - pd.Timedelta('1 hours')
t2 = edate + pd.Timedelta('1 hours')
file = cnst.MCS_15K# MCS_15K #_POINTS_DOM
msg = xr.open_dataarray(file)
try:
msg = msg.sel(time=slice(t1.strftime("%Y-%m-%dT%H"), t2.strftime("%Y-%m-%dT%H")))
except OverflowError:
return None
#print(prev_time.strftime("%Y-%m-%dT%H"), date.strftime("%Y-%m-%dT%H"))
pos = np.where((msg.values <= -40) ) #(msg.values >= 5) & (msg.values < 65)) # #
out = np.zeros_like(msg)
out[pos] = 1
out = np.sum(out, axis=0)
out[out>0]=1
# if np.sum(out>1) != 0:
# 'Stop!!!'
# pdb.set_trace()
msg = msg.sum(axis=0)*0
xout = msg.copy()
xout.name = 'probs'
xout.values = out
return xout
开发者ID:cornkle,项目名称:proj_CEH,代码行数:48,代码来源:composite_backtrack_ERA_map0.py
示例16: t_trend_slice
def t_trend_slice():
#file = '/users/global/cornkle/data/ERA-I monthly/ERA-WA-Monthly-2mTemp.nc'
file = '/localscratch/wllf030/cornkle/ERA-I/monthly/old/ERA-Int-Monthly-2mTemp.nc'
fpath = '/users/global/cornkle/figs/CLOVER/months/'
dam = xr.open_dataarray(file)
lower = 9
higher = 11
da = dam[(dam['time.month']>=lower) & (dam['time.month']<=higher)]
da = da.sel(longitude=slice(-18,51), latitude=slice(36, -37))
da = da.groupby('time.year').mean(axis=0)
lons = da.longitude
lats = np.flip(da.latitude.values, axis=0)
# define a function to compute a linear trend of a timeseries
def linear_trend(x):
#pf = np.polyfit(np.arange(len(x)), x, 1)
pf, slope, int, p, ind = mk.test(np.arange(len(x)),x.squeeze().values, eps=0.001, alpha=0.01, Ha='upordown')
# we need to return a dataarray or else xarray's groupby won't be happy
if ind == 1:
issig = slope
else:
issig = np.nan
return xr.DataArray(issig, )
# stack lat and lon into a single dimension called allpoints
stacked = da.stack(allpoints=['latitude','longitude'])
# apply the function over allpoints to calculate the trend at each point
trend = stacked.groupby('allpoints').apply(linear_trend)
# unstack back to lat lon coordinates
trend_unstacked = trend.unstack('allpoints')
trend_unstacked = trend_unstacked*10. # warming over decade
da2 = xr.DataArray(trend_unstacked, coords=[lats, lons], dims=['latitude', 'longitude'])
fp = fpath + 'ttrend_'+str(lower).zfill(2)+'-'+str(higher).zfill(2)+'.png'
up.quick_map_salem(da2, vmin=-0.4, vmax=0.4, cmap='RdBu_r', save=fp) #
plt.close('all')
开发者ID:cornkle,项目名称:proj_CEH,代码行数:48,代码来源:era_tgrad.py
示例17: create_ancils
def create_ancils():
dummy = xr.open_dataarray(dummy_grid)
ds = xr.Dataset(attrs=dummy.attrs)
dummy = dummy.isel(grid_longitude_t=slice(box[0], box[1]), grid_latitude_t=slice(box[2], box[3]))
files = glob.glob(ancils+'*.nc')
files
for f in files:
varsdat = xr.open_dataset(f, decode_times=False)
if 'pseudo' in varsdat.keys():
varsdat = varsdat.isel(rlon=slice(box[0], box[1]), rlat=slice(box[2], box[3]))
data = varsdat['field1391'].values[0, 0, :,:].squeeze() # time, plant type, y, x
if 'past' in f:
var = 'veg_past'
elif 'current' in f:
var = 'veg_current'
else:
print('Ancils not found')
return
ds[var] = xr.DataArray(data, coords={'false_latitude': dummy.grid_latitude_t.values,
'false_longitude': dummy.grid_longitude_t.values,
'true_latitude': (
['false_latitude', 'false_longitude'], dummy.latitude_t.values),
'true_longitude': (
['false_latitude', 'false_longitude'], dummy.longitude_t.values)},
dims=['false_latitude', 'false_longitude'])
if 'ht' in varsdat.keys():
varsdat = varsdat.isel(rlon=slice(box[0], box[1]), rlat=slice(box[2], box[3]))
data = varsdat['ht'].values[0, 0, :, :].squeeze() # time, plant type, y, x
ds['topo'] = xr.DataArray(data, coords={'false_latitude': dummy.grid_latitude_t.values,
'false_longitude': dummy.grid_longitude_t.values,
'true_latitude': (
['false_latitude', 'false_longitude'], dummy.latitude_t.values),
'true_longitude': (
['false_latitude', 'false_longitude'], dummy.longitude_t.values)},
dims=['false_latitude', 'false_longitude'])
ds.to_netcdf(out+'ancils/ancils_vera.nc')
开发者ID:cornkle,项目名称:proj_CEH,代码行数:48,代码来源:JASMIN_extract_script_VERA.py
示例18: timeline_trend_count
def timeline_trend_count():
msg_folder = cnst.GRIDSAT
fname = 'aggs/gridsat_WA_-70_monthly_count_-40base_1000km2.nc'
da = xr.open_dataarray(msg_folder + fname)
da = da.sel(lat=slice(4.5,8), lon=slice(-10, 15))
#da=da.sel(lat=slice(5,10))
#da[da==0]=np.nan
mean = da.mean(dim=['lat', 'lon'])
#mean = mean[(mean['time.month']==8)]
f= plt.figure(figsize=(10,6))
for i in range(3,6):
bla = mean[(mean['time.month'] == i)]
bla.plot(label=str(i), marker='o')
plt.title('Average number of pixels <= -70C, 4.5-8N')
plt.legend()
开发者ID:cornkle,项目名称:proj_CEH,代码行数:16,代码来源:gridsat_baseplots.py
示例19: timeline_trend_mean
def timeline_trend_mean():
msg_folder = '/users/global/cornkle/data/OBS/gridsat/gridsat_netcdf/'
fname = 'gridsat_WA_-70_monthly.nc'
da = xr.open_dataarray(msg_folder + fname)
da=da.sel(lat=slice(5,7), lon=slice(-17,20))
da[da==0]=np.nan
mean = da.mean(dim=['lat', 'lon'])
#mean = mean[(mean['time.month']==8)]
f= plt.figure(figsize=(10,6))
for i in range(4,6):
bla = mean[(mean['time.month'] == i)]
bla.plot(label=str(i), marker='o')
plt.title('Monthly mean temperature of pixels <= -40C, 11-18N')
plt.legend()
plt.ylim(-78,-71)
开发者ID:cornkle,项目名称:proj_CEH,代码行数:16,代码来源:cp4_baseplots.py
示例20: run
def run(fieldnum):
unmask_cube = fits.getdata(
'/data6/piyanat/projects/fg1p/mc_cubes_heraxx/hera331/'
'hera331_mc_cube_p{:03d}.fits'.format(fieldnum)
)
input_dir = '/data6/piyanat/projects/fg1p/masked_cubes_heraxx/m1/' \
'{:d}MHz/p{:03d}'.format(args.bw, fieldnum)
s_arr0 = np.empty((4, nbins, 2, 50))
s_arr1 = np.empty((4, nbins, 2, 50))
s_f0 = np.empty(nbins)
for binnum in range(nbins):
data_da = xr.open_dataarray(
'{:s}/masked_cube_hera331_p{:03d}_m1_bw{:d}MHz_bin{:02d}_fpad.nc'
.format(input_dir, fieldnum, args.bw, binnum)
)
fov_window = get_fov_window(data_da)
f_window = get_f_window(data_da)
full_mask = fov_window[None, :, :] * f_window[:, None, None]
f_mask = np.ones_like(full_mask, dtype=bool) * f_window[:, None, None]
fov_mask = np.ones_like(full_mask, dtype=bool) * fov_window[None, :, :]
s_arr1[0, binnum] = cal_pdf(data_da.values.ravel())
s_arr1[1, binnum] = cal_pdf(data_da.values[f_mask].ravel())
s_arr1[2, binnum] = cal_pdf(data_da.values[fov_mask].ravel())
s_arr1[3, binnum] = cal_pdf(data_da.values[full_mask].ravel())
ch_cut = slice(0 + (binnum * nf), (nf * 2) + (binnum * nf))
unmask_data = unmask_cube[::-1][ch_cut][::-1]
s_arr0[0, binnum] = cal_pdf(unmask_data.ravel())
s_arr0[1, binnum] = cal_pdf(unmask_data[f_mask].ravel())
s_arr0[2, binnum] = cal_pdf(unmask_data[fov_mask].ravel())
s_arr0[3, binnum] = cal_pdf(unmask_data[full_mask].ravel())
s_f0[binnum] = data_da.f0.values
s_ds = xr.Dataset(
{'original': (['cut', 'f0', 'val', 'val_bin'], s_arr0),
'fg_masked': (['cut', 'f0', 'val', 'val_bin'], s_arr1)},
coords={'cut': np.array(['none', 'freq', 'fov', 'all']),
'f0': s_f0,
'val': ['pdf', 'bin_center'],
'val_bin': np.arange(50)}
)
s_ds.to_netcdf('/data6/piyanat/projects/fg1p/stats/'
'{:d}MHz/pdf_hera331_masked_cube_p{:03d}_bw{:d}MHz.nc'
.format(args.bw, fieldnum, args.bw))
开发者ID:piyanatk,项目名称:sim,代码行数:47,代码来源:calculate_pdf_from_signal_cube.py
注:本文中的xarray.open_dataarray函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论