<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Blackout gap detection example notebook
This notebook shows
1. the usage of a blackout gap detection on a dataset of 500 CMLs and a 10 day period
2. two example CMLs with which have a very high and a very low number of blackout gaps
Step1:  CMLs with more than 35 minutes of blackout within the 10-day example data
Step2:  Investiagte blackouts for two CMLs over three years
Data is available from  this zenodo repository and has to be downloaded for the following analysis
Step3:  The data consists of RSL and TSL from two CMLs with one minutes resolution over three years and respective path-averaged reference data from RADKLIM-YW, a gauge-adjusted, climatologically correct weather radar product from the German Weather Service.
<ASSISTANT_TASK:>
Python Code:
import matplotlib.pyplot as plt
import numpy as np
import pycomlink as pycml
import xarray as xr
from tqdm import tqdm
import urllib.request
import io
import pycomlink.processing.blackout_gap_detection as blackout_detection
# Do show xarray.Dataset representation as text because gitlab/github
# do not (yet) render the html output correctly (or sometimes not at all...)
xr.set_options(display_style="text");
# get data from 500 CMLs and eleven days
data_path = pycml.io.examples.get_example_data_path()
cmls = xr.open_dataset(data_path + "/example_cml_data.nc")
# get recieved signal levels values without default fill values (-99 and -99.9)
rsl = cmls.rsl.where(cmls.rsl > -99).isel(channel_id=0)
gap_start_list = []
gap_end_list = []
gap_mask_list = []
for cml_id in tqdm(rsl.cml_id):
    # select individual RSL time series
    rsl_i = rsl.sel(cml_id=cml_id)
    # check whether RSL before and after gaps are below rsl_threshold
    # here we take -65 dBm as threshold as all CMLs have a median RSL > -60 dBm
    gap_start, gap_end = blackout_detection.get_blackout_start_and_end(
        rsl=rsl_i.values, rsl_threshold=-65
    )
    
    # create a mask for all gaps fullflilling the criteria above
    mask = blackout_detection.created_blackout_gap_mask_from_start_end_markers(
        rsl_i.values, gap_start, gap_end
    )
    mask_reverse = blackout_detection.created_blackout_gap_mask_from_start_end_markers(
        rsl_i.values[::-1], gap_end[::-1], gap_start[::-1]
    )
    mask = mask | mask_reverse[::-1]
    gap_start_list.append(gap_start)
    gap_end_list.append(gap_end)
    gap_mask_list.append(mask)
# parse gap starts, ends and the mask to one xarray dataset with the CML data
rsl = rsl.to_dataset()
rsl["gap_start"] = (["cml_id", "time"], gap_start_list)
rsl["gap_end_list"] = (["cml_id", "time"], gap_end_list)
rsl["mask"] = (["cml_id", "time"], gap_mask_list)
# get the CMLs with the many blackout minutes (here more than 35)
rsl_blackouts = rsl.isel(cml_id=(np.array(gap_mask_list).sum(axis=1) > 35))
rsl_blackouts
Explanation: Blackout gap detection example notebook
This notebook shows
1. the usage of a blackout gap detection on a dataset of 500 CMLs and a 10 day period
2. two example CMLs with which have a very high and a very low number of blackout gaps
End of explanation
# plot the three CMLs with most blackouts minutes (more than 35 minutes)
for cml_id in rsl_blackouts.cml_id:
    rsl_blackouts.rsl.sel(cml_id=cml_id, time="2018-05-13").plot(
        figsize=(10, 4), label="rsl"
    )
    (rsl_blackouts.mask.sel(cml_id=cml_id, time="2018-05-13") * -20).plot(
        label="detected blackout"
    )
    plt.legend()
Explanation: CMLs with more than 35 minutes of blackout within the 10-day example data
End of explanation
# load data e.g. with. curl -0` or `wget` 
!curl -O https://zenodo.org/record/6337557/files/blackout_example_cmls.nc
!curl -O https://zenodo.org/record/6337557/files/blackout_example_radar_reference.nc
# open dataset with xarray
cmls = xr.open_dataset('blackout_example_cmls.nc')
reference = xr.open_dataset('blackout_example_radar_reference.nc')
Explanation: Investiagte blackouts for two CMLs over three years
Data is available from  this zenodo repository and has to be downloaded for the following analysis
End of explanation
cmls
reference
# remove default rsl values
cmls["rsl"] = cmls.rsl.where(cmls.rsl > -99).isel(channel_id=0)
cmls["tsl"] = cmls.tsl.where(cmls.rsl > -99).isel(channel_id=0)
# define a plotting function
def plt_ts(rsl, tsl, ref, mask, start, ax=None):
    if ax is None:
        fig, ax = plt.subplots()
    ay = ax.twinx()
    ay.bar(np.arange(len(ref)), ref * 60, color="#045a8d", alpha=0.6)
    ay.set_ylim(-10, 180)
    ay.set_yticks([0, 50, 100, 150])
    ay.set_ylabel("rainfall intensity \n[mm/h]", color="#045a8d")
    ax.plot(tsl, color="#238b45", label="transmitted \nsignal level", lw=2)
    ax.plot(rsl, color="#cc4c02", label="recieved \nsignal level", lw=2)
    ax.plot((mask * -5) + 33.5, lw=4, color="black")
    ax.plot((mask + 100), color="black", lw=3, label="detected \nblackout gap")
    ax.set_ylim(-92, 30)
    ax.set_title("")
    ax.set_xticks([0, 30, 60, 90, 120, 150, 180, 210, 240])
    ax.set_xlabel("time [minutes]")
    ax.legend(
        loc="center left", ncol=3, 
    )
    ax.set_title(str(start.values)[0:10] + " " + str(start.values)[11:16])
    ax.set_yticks([20, 0, -20, -40, -60, -80])
    ax.set_ylabel("signal level\n[dBm]")
    plt.show()
for cml_id in cmls.cml_id:
    rsl = cmls.sel(cml_id=cml_id).rsl
    tsl = cmls.sel(cml_id=cml_id).tsl
    ref = reference.sel(cml_id=cml_id)
    
    # using the blackout gap detection as in the example above
    gap_start, gap_end = blackout_detection.get_blackout_start_and_end(
        rsl=rsl.values, rsl_threshold=-65
    )
    mask = blackout_detection.created_blackout_gap_mask_from_start_end_markers(
        rsl.values, gap_start, gap_end
    )
    mask_reverse = blackout_detection.created_blackout_gap_mask_from_start_end_markers(
        rsl.values[::-1], gap_end[::-1], gap_start[::-1]
    )
    mask = mask | mask_reverse[::-1]
    rsl = rsl.to_dataset()
    rsl["mask"] = ("time", mask)
    print(
        "For CML with the id "
        + str(cml_id.values)
        + " there are "
        + str(mask.sum())
        + " detected blackout minutes."
    )
    print("Plotting all detected gaps:")
    # remove gap_start times which are less than 60 minutes after another gap_start
    very_close_gaps = (
        np.diff(rsl.sel(time=gap_start).time.values) / 60000000000
    ).astype(int) < 60
    gap_time = rsl.sel(time=gap_start).time[~np.append(very_close_gaps, np.array(True))]
    # plot each gap
    for gap_start_time in gap_time:
        start = gap_start_time - np.timedelta64(2, "h")
        end = gap_start_time + np.timedelta64(2, "h")
        plt_ts(
            rsl=rsl.rsl.sel(time=slice(start, end)),
            tsl=tsl.sel(time=slice(start, end)),
            ref=ref.rainfall_amount.sel(time=slice(start, end)),
            mask=rsl.mask.sel(time=slice(start, end)),
            start=start,
        )
    print("#################################################\n")
Explanation: The data consists of RSL and TSL from two CMLs with one minutes resolution over three years and respective path-averaged reference data from RADKLIM-YW, a gauge-adjusted, climatologically correct weather radar product from the German Weather Service.
End of explanation
<END_TASK>