Ger radar rainfall along CML paths

import pycomlink as pycml
import xarray as xr
import matplotlib.pyplot as plt
import poligrain as plg

Load example data

data_path = pycml.io.examples.get_example_data_path()
cmls = xr.open_dataset(data_path + '/example_cml_data.nc')

ds_radar = xr.open_dataset(data_path + '/example_areal_reference_data.nc')

Calculate intersection weights

Note that the intersection weights are stored as sparse.arrays in a xarray.DataArray because the matrix of intersection weights for each CMLs, based on the radar grid, contains mostly zeros. Hence, we can save a lot of space. For large CML networks this is crucial because storing thousands of intersection weight matrices, one for each CML, easily eats up 10s of GBs of memory.

Note that this calculation is fairly fast, i.e. approx. 2 seconds for 500 CMLs.

%%time

da_intersect_weights = plg.spatial.calc_sparse_intersect_weights_for_several_cmls(
    x1_line=cmls.site_a_longitude.values,
    y1_line=cmls.site_a_latitude.values,
    x2_line=cmls.site_b_longitude.values,
    y2_line=cmls.site_b_latitude.values,
    cml_id=cmls.cml_id.values,
    x_grid=ds_radar.longitudes.values,
    y_grid=ds_radar.latitudes.values,
    grid_point_location='lower_left',
)

CPU times: user 2.39 s, sys: 29.4 ms, total: 2.42 s
Wall time: 2.44 s

da_intersect_weights

<xarray.DataArray (cml_id: 500, y: 190, x: 228)> Size: 155kB
<COO: shape=(500, 190, 228), dtype=float64, nnz=4842, fill_value=0.0>
Coordinates:
    x_grid   (y, x) float64 347kB 1.23 1.24 1.25 1.27 ... 4.14 4.15 4.17 4.18
    y_grid   (y, x) float64 347kB 56.88 56.88 56.88 56.88 ... 58.48 58.48 58.47
  * cml_id   (cml_id) <U3 6kB '0' '1' '2' '3' '4' ... '496' '497' '498' '499'
Dimensions without coordinates: y, x

xarray.DataArray

• cml_id: 500
• y: 190
• x: 228

• <COO: nnz=4842, fill_value=0.0>

  Format	coo
  Data Type	float64
  Shape	(500, 190, 228)
  nnz	4842
  Density	0.00022354570637119113
  Read-only	True
  Size	151.3K
  Storage ratio	0.00

• Coordinates: (3)
  • x_grid
    (y, x)
    float64
    1.23 1.24 1.25 ... 4.15 4.17 4.18

    array([[1.23, 1.24, 1.25, ..., 4.09, 4.1 , 4.11],
           [1.23, 1.24, 1.25, ..., 4.09, 4.1 , 4.11],
           [1.23, 1.24, 1.25, ..., 4.09, 4.1 , 4.11],
           ...,
           [1.17, 1.18, 1.2 , ..., 4.15, 4.17, 4.18],
           [1.17, 1.18, 1.2 , ..., 4.15, 4.17, 4.18],
           [1.17, 1.18, 1.2 , ..., 4.15, 4.17, 4.18]])

  • y_grid
    (y, x)
    float64
    56.88 56.88 56.88 ... 58.48 58.47

    array([[56.88, 56.88, 56.88, ..., 56.87, 56.87, 56.87],
           [56.89, 56.89, 56.89, ..., 56.88, 56.88, 56.88],
           [56.89, 56.89, 56.9 , ..., 56.89, 56.89, 56.89],
           ...,
           [58.46, 58.46, 58.46, ..., 58.46, 58.46, 58.46],
           [58.47, 58.47, 58.47, ..., 58.47, 58.47, 58.47],
           [58.48, 58.48, 58.48, ..., 58.48, 58.48, 58.47]])

  • cml_id
    (cml_id)
    <U3
    '0' '1' '2' ... '497' '498' '499'

    array(['0', '1', '2', ..., '497', '498', '499'], dtype='<U3')

• Indexes: (1)
  • cml_id
    PandasIndex

    PandasIndex(Index(['0', '1', '2', '3', '4', '5', '6', '7', '8', '9',
           ...
           '490', '491', '492', '493', '494', '495', '496', '497', '498', '499'],
          dtype='object', name='cml_id', length=500))

• Attributes: (0)

Plot intersection weights and CML path

Note that the lon-lat of the radar grid have been rounded in this example dataset and hence there are some strange pixel size, stemming from the non-equidistant lon-lat grids of the original radar data.

i = 22

fig, ax = plt.subplots(figsize=(5, 5))
ax.pcolormesh(
    ds_radar.longitudes.values,
    ds_radar.latitudes.values,
    da_intersect_weights.isel(cml_id=i).to_numpy()[:-1, :-1],
    shading='flat',
    cmap='Greys',
)

cml = cmls.isel(cml_id=i)
plt.ylim(
    min(cml.site_a_latitude.values, cml.site_b_latitude.values) - 0.1,
    max(cml.site_a_latitude.values, cml.site_b_latitude.values) + 0.1,
)
plt.xlim(
    min(cml.site_a_longitude.values, cml.site_b_longitude.values) - 0.1,
    max(cml.site_a_longitude.values, cml.site_b_longitude.values) + 0.1,
)

def plot_cml_lines(ds_cmls, ax):
    ax.plot(
        [ds_cmls.site_a_longitude, ds_cmls.site_b_longitude],
        [ds_cmls.site_a_latitude, ds_cmls.site_b_latitude],
        'r',
        linewidth=1,
    )

plot_cml_lines(ds_cmls=cmls.isel(cml_id=i), ax=plt.gca())

Get time series of radar rainfall along CML for each CML

This is internally done via a sparase.tensordot product, which is very fast, given that the intersection weights are also stored as sparse.arrays. The convenience funtion from pycomlink used here, also takes care of adding the correct coordinates to the xarray.DataArray that is returned.

Note that this calculation is very fast and can be done with a large number of CMLs and a much longer period of radar data.

%%time

da_radar_along_cmls = plg.spatial.get_grid_time_series_at_intersections(
    grid_data=ds_radar.rainfall_amount,
    intersect_weights=da_intersect_weights,
)

CPU times: user 2.41 s, sys: 1.29 s, total: 3.7 s
Wall time: 4.15 s

da_radar_along_cmls

<xarray.DataArray (time: 3168, cml_id: 500)> Size: 13MB
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])
Coordinates:
  * time     (time) datetime64[ns] 25kB 2018-05-10 ... 2018-05-20T23:55:00
  * cml_id   (cml_id) <U3 6kB '0' '1' '2' '3' '4' ... '496' '497' '498' '499'

xarray.DataArray

• time: 3168
• cml_id: 500

• 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

  array([[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]])

• Coordinates: (2)
  • time
    (time)
    datetime64[ns]
    2018-05-10 ... 2018-05-20T23:55:00

    array(['2018-05-10T00:00:00.000000000', '2018-05-10T00:05:00.000000000',
           '2018-05-10T00:10:00.000000000', ..., '2018-05-20T23:45:00.000000000',
           '2018-05-20T23:50:00.000000000', '2018-05-20T23:55:00.000000000'],
          dtype='datetime64[ns]')

  • cml_id
    (cml_id)
    <U3
    '0' '1' '2' ... '497' '498' '499'

    array(['0', '1', '2', ..., '497', '498', '499'], dtype='<U3')

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2018-05-10 00:00:00', '2018-05-10 00:05:00',
                   '2018-05-10 00:10:00', '2018-05-10 00:15:00',
                   '2018-05-10 00:20:00', '2018-05-10 00:25:00',
                   '2018-05-10 00:30:00', '2018-05-10 00:35:00',
                   '2018-05-10 00:40:00', '2018-05-10 00:45:00',
                   ...
                   '2018-05-20 23:10:00', '2018-05-20 23:15:00',
                   '2018-05-20 23:20:00', '2018-05-20 23:25:00',
                   '2018-05-20 23:30:00', '2018-05-20 23:35:00',
                   '2018-05-20 23:40:00', '2018-05-20 23:45:00',
                   '2018-05-20 23:50:00', '2018-05-20 23:55:00'],
                  dtype='datetime64[ns]', name='time', length=3168, freq=None))

  • cml_id
    PandasIndex

    PandasIndex(Index(['0', '1', '2', '3', '4', '5', '6', '7', '8', '9',
           ...
           '490', '491', '492', '493', '494', '495', '496', '497', '498', '499'],
          dtype='object', name='cml_id', length=500))

• Attributes: (0)

# The following show how you write the data to a NetCDF to use it somewhere else.
# This is how the path-averaged reference data in the example data directory was produced.
#
# da_radar_along_cmls.to_dataset(name='rainfall_amount').to_netcdf(
#    '../pycomlink/io/example_data/example_path_averaged_reference_data.nc',
#    encoding={'rainfall_amount': {'zlib': True}}
#)

Plot radar along CML vs. TRSL

We set non-meaningfull RSL and TSL values to NaN.

Here we only plot radar rain rates vs TRSL. For plotting CML-derived rain rates, processing has to be done first, which is e.g. shown in the notebook with the basic processing workflow. But also the comparison with TRSL already shows how good CML-data corresponds to rainfall data.

(we might add this processing here later…)

cmls['tsl'] = cmls.tsl.where(cmls.tsl < 100)
cmls['rsl'] = cmls.rsl.where(cmls.rsl > -99)

cmls['trsl'] = cmls.tsl - cmls.rsl

for i in range(10):
    fig, axs = plt.subplots(2, 1, figsize=(14,4), sharex=True)
    da_radar_along_cmls.isel(cml_id=i).plot(ax=axs[0])
    cmls.isel(cml_id=i).trsl.plot.line(x='time', ax=axs[1])
    axs[0].set_xlabel('')
    axs[0].set_ylabel('radar 5-min rainfall sum\nalong CML\n(mm/h)')
    axs[1].set_title('')
    axs[1].set_ylabel('trsl\n(dB)')