Python xarray.open_dataarray函数代码示例

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()

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()))

    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)

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

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')

 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)

 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)

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)

    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

 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)

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')

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)

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')

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()

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