Merge multisensor data

This module groups all the functions used to merge together different data sources such as CML, Gauges or Radar.

Check the Sample data doc page to see how cml and gauges datasets are structured and what variables they contain.

cml, gauges = open_cml_sample(),open_gauge_sample()

Assign nearest gauge using link center approximation

To assign a gauge to each link, we approximate the link’s position by its center point (midpoint of the two antenna sites), then find the closest gauge to that center.

The distance is computed using the haversine formula, which gives the great-circle distance between two points \((\phi_1, \lambda_1)\) and \((\phi_2, \lambda_2)\) on a sphere:

\[d = 2R \arcsin\left(\sqrt{\sin^2\frac{\Delta\phi}{2} + \cos\phi_1 \cos\phi_2 \sin^2\frac{\Delta\lambda}{2}}\right)\]

where \(\Delta\phi = \phi_2 - \phi_1\), \(\Delta\lambda = \lambda_2 - \lambda_1\) (in radians), and \(R \approx 6371\,\text{km}\) is the mean Earth radius.

This is a spherical approximation — it assumes the Earth is a perfect sphere rather than using an ellipsoidal model (e.g. WGS-84). The resulting error is at most ~0.5%, which is negligible here given that using the link center is already an approximation of the true link path.

# Convert to radians — shapes: (n_links,1) and (1,n_gauges)
lat1 = np.radians(links_df["lat"].values[:, None])
lon1 = np.radians(links_df["lon"].values[:, None])
lat2 = np.radians(gauge_coords["lat"].values[None, :])
lon2 = np.radians(gauge_coords["lon"].values[None, :])

print(lat1.shape, lat2.shape, lat2[0] - lat1[0])

(126, 1) (1, 8) [0.00125581 0.00125985 0.00164869 0.00163818 0.0008165  0.00126027
 0.00117194 0.00109306]

dlat = lat2 - lat1
dlon = lon2 - lon1
print(dlat.shape, dlon.shape)
dlat[0]

(126, 8) (126, 8)

array([0.00125581, 0.00125985, 0.00164869, 0.00163818, 0.0008165 ,
       0.00126027, 0.00117194, 0.00109306])

a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
dist_km = 6371 * 2 * np.arcsin(np.sqrt(a))  # (n_links, n_gauges)

sd = distances.sort_values().to_frame(name="Geodesic").copy()
sd["Haversine"] = sorted(dist_km[0])
sd["Diff (km)"] = (sd["Geodesic"] - sd["Haversine"]).abs()
sd

	Geodesic	Haversine	Diff (km)
id
Lycee_NYALLA	5.612400	5.637986	0.025587
Ecole_Saint_Andre_PK11	8.029850	8.073887	0.044037
SOCATUR	8.475717	8.504042	0.028325
Meteo_IUT	9.381050	9.415967	0.034917
Campus_2_Universite	9.526730	9.560472	0.033742
Piscine_CITE_SIC	9.896087	9.922225	0.026138
Eglise_Pie_X_PK14	10.461697	10.519682	0.057985
Hopital_des_Soeurs_Logpom	10.654718	10.710358	0.055640

# Find nearest gauge per link
idx = np.argmin(dist_km, axis=1)
min_dist = dist_km[np.arange(len(idx)), idx]
gauge_coords.index[idx[0]], min_dist[0]

('Lycee_NYALLA', np.float64(5.637986481941067))

nearest = pd.DataFrame({
    "id": gauge_coords.index[idx],
    "distance_km": min_dist,
}, index=links_df.index)
nearest.head(2)

	id	distance_km
cml_id
3.984686N-9.789517E	Lycee_NYALLA	5.637986
3.985691N-9.801648E	Lycee_NYALLA	6.188639

gauge_n_dist.head()

CPU times: user 1.17 ms, sys: 0 ns, total: 1.17 ms
Wall time: 975 μs

	id	distance_km
cml_id
3.984686N-9.789517E	Lycee_NYALLA	5.637986
3.985691N-9.801648E	Lycee_NYALLA	6.188639
3.996996N-9.761328E	Lycee_NYALLA	3.949761
4.002892N-9.746861E	Lycee_NYALLA	4.078766
4.007006N-9.767519E	Lycee_NYALLA	2.733009

source

assign_nearest_gauge_to_link_center


def assign_nearest_gauge_to_link_center(
    cml:Dataset, # CML dataset containing link information
    gauges:Dataset, # Dataset containing gauges coordinates
    gauge_id:str='id', # Name of the dimension of the gauge identifier in the gauges dataset
    gauge_lat:str='lat', # Name of the latitude coordinate in the gauges dataset
    gauge_lon:str='lon', # Name of the longitude coordinate in the gauges dataset
    metadata_only:bool=False, # Assign only nearest gauge metadata
)->Dataset: # CML dataset merged with nearest gauge data and metadata

Assign the nearest gauge to each CML link based on geodesic distance from gauge to link center coordinates.

cml_with_gauges = assign_nearest_gauge_to_link_center(cml, gauges)
cml_with_gauges

cml_with_gauges_metadata = assign_nearest_gauge_to_link_center(cml, gauges, metadata_only=True)
cml_with_gauges

You can create an interactive version of this map using our interactive map generation tools.

Sample data generation

As we will need sample data from different processing steps to help with development and illustration, we will save the generated wet classification as a sample.

cml_with_gauges = assign_nearest_gauge_to_link_center(cml, gauges)
cml_with_gauges.to_netcdf(get_cml_sample_fp("with_nearest_gauge"), encoding={v: {"zlib": True, "complevel": 5} for v in cml_with_gauges.data_vars.keys()})