|
| 1 | +# conda activate venv_kde_xarray |
| 2 | + |
| 3 | +## Written on 18/01/2021 by Laura Gomez Navarro |
| 4 | + |
| 5 | +# Modules: |
| 6 | + |
| 7 | +from glob import glob |
| 8 | +import xarray as xr |
| 9 | +import numpy as np |
| 10 | +from scipy.stats import gaussian_kde |
| 11 | + |
| 12 | +def rem_nans(ds): |
| 13 | + """ |
| 14 | + This renders lon and lat variables without nans for the last timestep. |
| 15 | + """ |
| 16 | + bad_indices = np.isnan(ds['lon'][:,-1]) | np.isnan(ds['lat'][:,-1]) |
| 17 | + good_indices = ~bad_indices |
| 18 | + lon_end_nonans = ds['lon'][:,-1][good_indices] |
| 19 | + lat_end_nonans = ds['lat'][:,-1][good_indices] |
| 20 | + |
| 21 | + return lon_end_nonans, lat_end_nonans |
| 22 | + |
| 23 | +def kde_vals(lon_end_nonans, lat_end_nonans): |
| 24 | + """ |
| 25 | + TO DO |
| 26 | + """ |
| 27 | + x = lon_end_nonans.copy() |
| 28 | + y = lat_end_nonans.copy() |
| 29 | + |
| 30 | + xy = np.vstack([x, y]) |
| 31 | + z = gaussian_kde(xy)(xy) |
| 32 | + |
| 33 | + idx = z.argsort() |
| 34 | + x, y, z = x[idx], y[idx], z[idx] |
| 35 | + |
| 36 | + return x, y, z |
| 37 | + |
| 38 | +def kde_parcels(nfile): |
| 39 | + """ |
| 40 | + Inputs |
| 41 | + nfile : filedir + filename |
| 42 | +
|
| 43 | + Outputs |
| 44 | + KDE for end fields of lon and lat |
| 45 | + """ |
| 46 | + |
| 47 | + ds = xr.open_dataset(nfile) |
| 48 | + |
| 49 | + # Remove nans: |
| 50 | + lon_end_nonans, lat_end_nonans = rem_nans(ds) |
| 51 | + |
| 52 | + kde_x, kde_y, kde_z = kde_vals(lon_end_nonans, lat_end_nonans) |
| 53 | + |
| 54 | + savedir = '/data/oceanparcels/output_data/data_LauraGN/outputs_parcels/kde_calcs/' |
| 55 | + savename = savedir + 'KDE_' + nfile.split('/')[-1].split('.nc')[0] + '.npz' |
| 56 | + np.savez(savename, kde_x=kde_x, kde_y=kde_y, kde_z=kde_z) |
| 57 | + |
| 58 | + return kde_x, kde_y, kde_z |
| 59 | + |
| 60 | +def kde_parcels_bws(nfile, bw_set): |
| 61 | + """ |
| 62 | + bw_set = np.arange(0.01, 0.26, 0.01) |
| 63 | + """ |
| 64 | + ds = xr.open_dataset(nfile) |
| 65 | + |
| 66 | + # Remove nans: |
| 67 | + lon_end_nonans, lat_end_nonans = rem_nans(ds) |
| 68 | + |
| 69 | + x = lon_end_nonans.copy() |
| 70 | + y = lat_end_nonans.copy() |
| 71 | + |
| 72 | + xy = np.vstack([x, y]) |
| 73 | + |
| 74 | + KDE_dict = {} |
| 75 | + for ii in range(0, len(bw_set)): |
| 76 | + KDE_dict["z%02d" %ii] = gaussian_kde(xy, bw_method=bw_set[ii])(xy) |
| 77 | + |
| 78 | + # Saving: |
| 79 | + savedir = '/data/oceanparcels/output_data/data_LauraGN/outputs_parcels/kde_calcs/' |
| 80 | + savename = savedir + 'KDE_' + nfile.split('/')[-1].split('.nc')[0] + '_bws.npz' |
| 81 | + np.savez(savename, bw_set=bw_set, KDE_dict=KDE_dict) |
| 82 | + |
0 commit comments