Gaussian Kernel Density Estimation (method 2: with land mask)#
This code uses the KDE function from sci-kit learn package to find the smoothing curve over the 2D histogram of particle counts in the horizontal plane and plots it over a grid with the same resolution as the particles’ trajectories.
It takes into account the grid cells that are over land to avoid unrealistic results.
-Author: Jimena Medina Rubio
-Created on: 18/03/2023
0. Imports and package versions#
[1]:
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import cmocean
import cmocean.cm as cmo
import pandas as pd
from scipy.interpolate import griddata
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV
1. Loading the data#
[2]:
path="../../Simulations/toy_data_01.nc"
ds=xr.open_dataset(path)
print(ds.keys)
<bound method Mapping.keys of <xarray.Dataset> Size: 627kB
Dimensions: (traj: 144, obs: 121)
Dimensions without coordinates: traj, obs
Data variables:
trajectory (traj, obs) float64 139kB ...
time (traj, obs) datetime64[ns] 139kB ...
lat (traj, obs) float32 70kB ...
lon (traj, obs) float32 70kB ...
z (traj, obs) float32 70kB ...
U (traj, obs) float32 70kB ...
V (traj, obs) float32 70kB ...
Attributes:
feature_type: trajectory
Conventions: CF-1.6/CF-1.7
ncei_template_version: NCEI_NetCDF_Trajectory_Template_v2.0
parcels_version: 2.3.1.dev20+g92f2fb90
parcels_mesh: spherical>
2. Definition of functions#
[3]:
def land_mask_interpolated(land_mask, ds, bins_x, bins_y):
"""
Interpolates the input land_mask from a netCDF file to the desired grid on which the KDE function
is displayed. The number of grid cells in the zonal & meridional direction is bins_x & bins_y respectively
"""
#use boolean indexing to select the values of longitude & latitude inside the domain
land_mask = land_mask.sel(lat=slice(ds['lat'].min(), ds['lat'].max()),
lon=slice(ds['lon'].min(), ds['lon'].max()))
#define fine grid coordinates
fine_x, fine_y = np.meshgrid(np.arange(0, np.shape(land_mask)[1]), np.arange(0, np.shape(land_mask)[0]))
#define coarse grid coordinates with binx_x & bins_y cells in the x/y direction
coarse_x, coarse_y = np.meshgrid(np.linspace(0, np.shape(land_mask)[1], bins_x), np.linspace(0, np.shape(land_mask)[0], bins_y))
#flatten fine grid coordinates and mask
flat_x, flat_y, flat_mask = fine_x.flatten(), fine_y.flatten(), land_mask.stack(z=('lat', 'lon'))
#interpolate mask onto coarse grid
coarse_mask = griddata((flat_x, flat_y), flat_mask, (coarse_x, coarse_y), method='nearest')
#set NaN values in ocean
coarse_mask[coarse_mask==100]=np.nan
#create boolean mask in the same grid as used in kde with True in ocean
ocean=np.isnan(coarse_mask)
return ocean
[4]:
def kde_landmask(ds, bins_x, bins_y, ocean):
"""
Calculates the smoothing Gaussian function of the longitude and latitude samples from the particles' locations
at each observation from the Dataset 'ds'. The bandwidth of such function is fitted to the samples and the
resulting heat map is shown over an adjustable grid controlled by the parameters 'binx_x' & 'bins_y'.
The input land mask 'ocean' is a boolean mask that restricts the KDE values to the ocean.
"""
#initialize an empty list to store the KDE values for each observation
kde_values = []
#loop over each observation in the dataset
for i in range(ds.obs.size):
#extract the lon/lat values for the current observation
lon_lat = np.vstack([ds.lat[:,i], ds.lon[:,i]]).T
#create a 2D grid of lat/lon values to evaluate the KDE on based on
x, y = np.meshgrid(np.linspace(ds['lon'].min(), ds['lon'].max(), bins_x), np.linspace(ds['lat'].min(), ds['lat'].max(), bins_y))
#flatten grid & the ocean mask
xy = np.column_stack([y.ravel(), x.ravel()])
ocean=ocean.ravel()
#only consider points from xy grid in ocean
xy = xy[ocean]
#initialize a KDE object & fit it to the lon/lat values from the trajectories
kde = KernelDensity( kernel="gaussian", algorithm="ball_tree")
#suggest a range of possible values of bandwidth & find best fit
bandwidth = np.arange(0.02, 1, 0.05)
grid = GridSearchCV(kde, {'bandwidth': bandwidth})
#find KDE for the given lat/lon values
grid.fit(lon_lat)
kde = grid.best_estimator_
#evaluate the KDE on the 2D grid created & get log-likelihood
log_density = kde.score_samples(xy)
#obtain probability at each grid cell
density = np.exp(log_density)
#normalise results
density /= density.sum()
#initialise output matrix
z = np.zeros(x.shape)
#only keep values in the ocean
z.ravel()[ocean] = density
z.resize(x.shape)
# append the KDE values of each grid cell to the list
kde_values.append(z)
#combine the KDE values into a DataArray
kde_values = np.stack(kde_values)
kde_da = xr.DataArray(kde_values,
dims=("obs", "lat", "lon"),
coords={"obs": ds.obs.values,
"lat": y[:, 0],
"lon": x[0, :]},
name="%")
#print the used bandwidth
print("optimal bandwidth: " + "{:.2f}".format(kde.bandwidth))
#compute the cumulative sum of the particle distribution over time
kde_total=np.sum(kde_da, axis=0)
#obtain normalised result in percentage
kde_total = kde_total*100/np.nansum(kde_total)
#set to NaN zero value
kde_total=kde_total.where(kde_total!= 0, np.nan)
return kde_total
[5]:
def probability_map(probability, xlim, ylim, title):
""" All-included plot of the desired domain specified by xlim & ylim """
fig=plt.figure(figsize=(13,6))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree(central_longitude=0.0))
#create grid with [min, max] values of lons & lats
ax.set_xlim(xlim)
ax.set_ylim(ylim)
#plot coastlines
ax.coastlines(resolution='10m')
ax.add_feature(cartopy.feature.LAND, facecolor='grey')
#draw grid lines
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gl.right_labels = False
gl.top_labels = False
#plot probability results
probability.plot(ax=ax, cmap=cmo.matter)
plt.title(title)
return plt.show()
3. Results#
3.1 Obtaining land mask to exclude from KDE analysis#
Method: imput fine resolution land mask and interpolate it to the grid used to display the KDE on. Download land sea mask from NASA website.
[6]:
#import land mask
land=xr.open_dataset('/Users/lienzo/Downloads/IMERG_land_sea_mask.nc')
print(land.keys())
land['landseamask'].plot()
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/file_manager.py:211, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
210 try:
--> 211 file = self._cache[self._key]
212 except KeyError:
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/lru_cache.py:56, in LRUCache.__getitem__(self, key)
55 with self._lock:
---> 56 value = self._cache[key]
57 self._cache.move_to_end(key)
KeyError: [<class 'netCDF4._netCDF4.Dataset'>, ('/Users/lienzo/Downloads/IMERG_land_sea_mask.nc',), 'r', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False)), 'bcff2111-a4b9-402c-a696-ac8b52719089']
During handling of the above exception, another exception occurred:
FileNotFoundError Traceback (most recent call last)
Cell In[6], line 2
1 #import land mask
----> 2 land=xr.open_dataset('/Users/lienzo/Downloads/IMERG_land_sea_mask.nc')
3 print(land.keys())
4 land['landseamask'].plot()
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/api.py:760, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, create_default_indexes, inline_array, chunked_array_type, from_array_kwargs, backend_kwargs, **kwargs)
748 decoders = _resolve_decoders_kwargs(
749 decode_cf,
750 open_backend_dataset_parameters=backend.open_dataset_parameters,
(...) 756 decode_coords=decode_coords,
757 )
759 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 760 backend_ds = backend.open_dataset(
761 filename_or_obj,
762 drop_variables=drop_variables,
763 **decoders,
764 **kwargs,
765 )
766 ds = _dataset_from_backend_dataset(
767 backend_ds,
768 filename_or_obj,
(...) 779 **kwargs,
780 )
781 return ds
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/netCDF4_.py:682, in NetCDF4BackendEntrypoint.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, format, clobber, diskless, persist, auto_complex, lock, autoclose)
660 def open_dataset(
661 self,
662 filename_or_obj: T_PathFileOrDataStore,
(...) 679 autoclose=False,
680 ) -> Dataset:
681 filename_or_obj = _normalize_path(filename_or_obj)
--> 682 store = NetCDF4DataStore.open(
683 filename_or_obj,
684 mode=mode,
685 format=format,
686 group=group,
687 clobber=clobber,
688 diskless=diskless,
689 persist=persist,
690 auto_complex=auto_complex,
691 lock=lock,
692 autoclose=autoclose,
693 )
695 store_entrypoint = StoreBackendEntrypoint()
696 with close_on_error(store):
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/netCDF4_.py:468, in NetCDF4DataStore.open(cls, filename, mode, format, group, clobber, diskless, persist, auto_complex, lock, lock_maker, autoclose)
464 kwargs["auto_complex"] = auto_complex
465 manager = CachingFileManager(
466 netCDF4.Dataset, filename, mode=mode, kwargs=kwargs
467 )
--> 468 return cls(manager, group=group, mode=mode, lock=lock, autoclose=autoclose)
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/netCDF4_.py:398, in NetCDF4DataStore.__init__(self, manager, group, mode, lock, autoclose)
396 self._group = group
397 self._mode = mode
--> 398 self.format = self.ds.data_model
399 self._filename = self.ds.filepath()
400 self.is_remote = is_remote_uri(self._filename)
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/netCDF4_.py:477, in NetCDF4DataStore.ds(self)
475 @property
476 def ds(self):
--> 477 return self._acquire()
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/netCDF4_.py:471, in NetCDF4DataStore._acquire(self, needs_lock)
470 def _acquire(self, needs_lock=True):
--> 471 with self._manager.acquire_context(needs_lock) as root:
472 ds = _nc4_require_group(root, self._group, self._mode)
473 return ds
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/contextlib.py:137, in _GeneratorContextManager.__enter__(self)
135 del self.args, self.kwds, self.func
136 try:
--> 137 return next(self.gen)
138 except StopIteration:
139 raise RuntimeError("generator didn't yield") from None
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/file_manager.py:199, in CachingFileManager.acquire_context(self, needs_lock)
196 @contextlib.contextmanager
197 def acquire_context(self, needs_lock=True):
198 """Context manager for acquiring a file."""
--> 199 file, cached = self._acquire_with_cache_info(needs_lock)
200 try:
201 yield file
File ~/miniforge3/envs/lagrangian_diags/lib/python3.11/site-packages/xarray/backends/file_manager.py:217, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
215 kwargs = kwargs.copy()
216 kwargs["mode"] = self._mode
--> 217 file = self._opener(*self._args, **kwargs)
218 if self._mode == "w":
219 # ensure file doesn't get overridden when opened again
220 self._mode = "a"
File src/netCDF4/_netCDF4.pyx:2521, in netCDF4._netCDF4.Dataset.__init__()
File src/netCDF4/_netCDF4.pyx:2158, in netCDF4._netCDF4._ensure_nc_success()
FileNotFoundError: [Errno 2] No such file or directory: '/Users/lienzo/Downloads/IMERG_land_sea_mask.nc'
[ ]:
#choose bin size of the output grid
bins_x=100
bins_y=100
#interpolate land mask to output grid
ocean=land_mask_interpolated(land['landseamask'], ds, bins_x, bins_y)
#calculate KDE over output grid
kde_results=kde_landmask(ds, bins_x, bins_y, ocean)
4. Plotting#
[ ]:
#choose the limits of the x & y axis of the graph
xlim=[10, 47]
ylim= [-47, -25]
#plot results
probability_map(kde_results, xlim, ylim, 'Kernel Distribution Estimation')
[ ]: