<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Table of Contents
<p><div class="lev1 toc-item"><a href="#Load-CML-example-data" data-toc-modified-id="Load-CML-example-data-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Load CML example data</a></div><div class="lev1 toc-item"><a href="#Do-a-simple-standard-processing-to-get-rain-rates-for-each-CML" data-toc-modified-id="Do-a-simple-standard-processing-to-get-rain-rates-for-each-CML-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Do a simple standard processing to get rain rates for each CML</a></div><div class="lev1 toc-item"><a href="#Do-IDW-interpolation-of-CML-rain-rates" data-toc-modified-id="Do-IDW-interpolation-of-CML-rain-rates-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>Do IDW interpolation of CML rain rates</a></div><div class="lev2 toc-item"><a href="#Initialize-interpolator" data-toc-modified-id="Initialize-interpolator-31"><span class="toc-item-num">3.1&nbsp;&nbsp;</span>Initialize interpolator</a></div><div class="lev2 toc-item"><a href="#Perform-interpolation-for-all-time-steps" data-toc-modified-id="Perform-interpolation-for-all-time-steps-32"><span class="toc-item-num">3.2&nbsp;&nbsp;</span>Perform interpolation for all time steps</a></div>
Step1:  Load CML example data
Coordinates mimic the real network topology but are fake
Step2:  Do a simple standard processing to get rain rates for each CML
Step3:  Do IDW interpolation of CML rain rates
The ComlinkGridInterpolator takes a PointsToGridInterpolator object as argument, which is used for the interpolation of each time step. You can pass config arguments to the initialization of the PointsToGridInterpolator.
Currently only the IDW interpolator IdWKdtreeInterpolator which subclasses PointsToGridInterpolator is available. A Kriging version is already implemented but does not work reliably.
Initialize interpolator
resolution is used to generate a grid using a bounding box aroudn all CMLs if no x- and y-grid are supplied.
Currently CML rain rates are averaged to hourly data before interpolating.
Step4:  Perform interpolation for all time steps
Step5:  Calculate CML coverage mask
Coverage for 0.05 degree coverage around CMLs.
Note
Step6:  Coverage for 0.1 degree coverage around CMLs.
Step7:  Plot CML rainfall sum and apply coverage map
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import pycomlink as pycml
import matplotlib.pyplot as plt
from tqdm import tqdm
Explanation: Table of Contents
<p><div class="lev1 toc-item"><a href="#Load-CML-example-data" data-toc-modified-id="Load-CML-example-data-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Load CML example data</a></div><div class="lev1 toc-item"><a href="#Do-a-simple-standard-processing-to-get-rain-rates-for-each-CML" data-toc-modified-id="Do-a-simple-standard-processing-to-get-rain-rates-for-each-CML-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Do a simple standard processing to get rain rates for each CML</a></div><div class="lev1 toc-item"><a href="#Do-IDW-interpolation-of-CML-rain-rates" data-toc-modified-id="Do-IDW-interpolation-of-CML-rain-rates-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>Do IDW interpolation of CML rain rates</a></div><div class="lev2 toc-item"><a href="#Initialize-interpolator" data-toc-modified-id="Initialize-interpolator-31"><span class="toc-item-num">3.1&nbsp;&nbsp;</span>Initialize interpolator</a></div><div class="lev2 toc-item"><a href="#Perform-interpolation-for-all-time-steps" data-toc-modified-id="Perform-interpolation-for-all-time-steps-32"><span class="toc-item-num">3.2&nbsp;&nbsp;</span>Perform interpolation for all time steps</a></div>
End of explanation
cml_list = pycml.io.examples.get_75_cmls()
fig, ax = plt.subplots()
for cml in cml_list:
    cml.plot_line(ax=ax, color='k')
Explanation: Load CML example data
Coordinates mimic the real network topology but are fake
End of explanation
for cml in tqdm(cml_list):
    window_length = 60
    threshold = 1.0
    cml.process.wet_dry.std_dev(window_length=window_length, threshold=threshold)
    cml.process.baseline.linear()
    cml.process.baseline.calc_A()
    cml.process.A_R.calc_R()
Explanation: Do a simple standard processing to get rain rates for each CML
End of explanation
cml_interp = pycml.spatial.interpolator.ComlinkGridInterpolator(
    cml_list=cml_list,
    resolution=0.01,
    interpolator=pycml.spatial.interpolator.IdwKdtreeInterpolator())
Explanation: Do IDW interpolation of CML rain rates
The ComlinkGridInterpolator takes a PointsToGridInterpolator object as argument, which is used for the interpolation of each time step. You can pass config arguments to the initialization of the PointsToGridInterpolator.
Currently only the IDW interpolator IdWKdtreeInterpolator which subclasses PointsToGridInterpolator is available. A Kriging version is already implemented but does not work reliably.
Initialize interpolator
resolution is used to generate a grid using a bounding box aroudn all CMLs if no x- and y-grid are supplied.
Currently CML rain rates are averaged to hourly data before interpolating.
End of explanation
ds = cml_interp.loop_over_time()
ds
fig, ax = plt.subplots(3, 3, sharex=True, sharey=True, figsize=(12,12))
for i, axi in enumerate(ax.flat):
    for cml in cml_list:
        cml.plot_line(ax=axi, color='k')
    
    pc = axi.pcolormesh(ds.lon,
                        ds.lat,
                        ds.R.isel(time=20+i), 
                        cmap=plt.get_cmap('BuPu', 8), 
                        vmin=0,
                        vmax=20)
    axi.set_title(cml_interp.df_cmls.index[20+i])
    
fig.subplots_adjust(right=0.9)
cbar_ax = fig.add_axes([0.95, 0.15, 0.02, 0.7])
fig.colorbar(pc, cax=cbar_ax, label='Hourly rainfall sum in mm');
Explanation: Perform interpolation for all time steps
End of explanation
cml_coverage_mask = pycml.spatial.coverage.calc_coverage_mask(
    cml_list=cml_list, 
    xgrid=ds.lon.values,
    ygrid=ds.lat.values,
    max_dist_from_cml=0.05)
fig, ax = plt.subplots()
for cml in cml_list:
        cml.plot_line(ax=ax, color='k')
ax.pcolormesh(ds.lon, ds.lat, cml_coverage_mask, cmap='gray');
Explanation: Calculate CML coverage mask
Coverage for 0.05 degree coverage around CMLs.
Note: Calculating coverage using lon-lat and degrees does result in distortions. In the future this will be done using a area preserving reprojection of the lon-lat coordinates before calculating coverage.
End of explanation
cml_coverage_mask = pycml.spatial.coverage.calc_coverage_mask(
    cml_list=cml_list, 
    xgrid=ds.lon.values,
    ygrid=ds.lat.values,
    max_dist_from_cml=0.1)
fig, ax = plt.subplots()
for cml in cml_list:
        cml.plot_line(ax=ax, color='k')
ax.pcolormesh(ds.lon, ds.lat, cml_coverage_mask, cmap='gray');
Explanation: Coverage for 0.1 degree coverage around CMLs.
End of explanation
fig, ax = plt.subplots()
for cml in cml_list:
    cml.plot_line(ax=ax, color='k')
pc = ax.pcolormesh(
    ds.lon, 
    ds.lat,
    ds.R.sum(dim='time').where(cml_coverage_mask),
    cmap=plt.get_cmap('BuPu', 32))
plt.colorbar(pc, label='rainfall sum in mm');
Explanation: Plot CML rainfall sum and apply coverage map
End of explanation
<END_TASK>