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

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

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

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

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)

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