Quickstart guide functional api

This example is designed as quickstart guide for the functional api of pylawr.

Introduction

In the following, we will use the functional api to process netCDF data from the HHG X-band weather radar. These steps are almost the same as for the low-level api.

The following processing steps are: 1. Read in data with pre-defined data handler for netcdf data 2. Remove noise 3. Filter out clutter 4. Interpolate missing values 5. Plot the data

import pylawr.functions.input as input_funcs
import pylawr.functions.transform as transform_funcs
import pylawr.functions.plot as plot_funcs

from pylawr.grid import PolarGrid

1. Read in data

To read in the ascii data, we will use a pre-defined function to read in HHG X-band weather radar netCDF data. This functions takes a file name and returns read reflectivity. For this, we need to open a file handler. We need to set the grid manually, because we cannot extract the grid from given netCDF data at the moment.

file_path = '../tests/data/lawr_l0_example.nc'
reflectivity, _ = input_funcs.read_lawr_nc_level0(file_path)
grid = PolarGrid()
reflectivity = reflectivity.lawr.set_grid_coordinates(grid)

print(reflectivity)

<xarray.DataArray 'dbz' (time: 1, azimuth: 360, range: 333)>
array([[[-25.083254,  43.722687,  35.157005, ...,  21.714653,
          21.636692,  21.650703],
        [-25.132263,  43.60427 ,  37.33281 , ...,  21.870094,
          22.084316,  21.731308],
        [-25.14114 ,  43.791496,  37.9972  , ...,  21.80893 ,
          21.812252,  21.956898],
        ...,
        [-25.39464 ,  43.42102 ,  37.75023 , ...,  21.459068,
          22.672329,  22.397419],
        [-25.292507,  43.703377,  36.075798, ...,  21.719315,
          22.868332,  22.19233 ],
        [-25.27496 ,  43.455482,  36.701538, ...,  21.591888,
          22.23751 ,  21.603336]]], dtype=float32)
Coordinates:
  * time     (time) datetime64[ns] 2016-06-07T16:21:30
  * azimuth  (azimuth) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
  * range    (range) float64 29.98 89.94 149.9 ... 1.982e+04 1.988e+04 1.994e+04
Attributes: (12/17)
    units:                dBZ
    standard_name:        equivalent_reflectivity_factor
    short_name:           dBZ
    Convections:          CF-1.0
    History:              2016-06-08 09:32 generated from orignal data
    Sample_Frequency_Hz:  5000000
    ...                   ...
    institution:          Meteorological Institute, Hamburg, Germany
    title:                LAWR data
    latitude:             53.56833
    longitude:            9.973997
    corr_offset:          29.5364
    tags:                 beam-expansion-corr

print(grid)

<pylawr.grid.polar.PolarGrid object at 0x7fca9251b2d0>

2. Remove noise

As next step, we want to remove existing noise from reflectivity. For this, we will use remove_noise from pylawr.functions.transform, which fits a NoiseRemover and uses this remover to subtract the noise.

reflectivity, noise_remover = transform_funcs.remove_noise(reflectivity)

print(reflectivity)
print('Determined noiselevel: {0:.3e}'.format(noise_remover.noiselevel))

<xarray.DataArray 'z' (time: 1, azimuth: 360, range: 333)>
array([[[2.75846010e-03, 2.35650633e+04, 3.27868306e+03, ...,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
        [2.72364866e-03, 2.29312039e+04, 5.41103340e+03, ...,
         3.61638714e+00, 1.04837054e+01, 0.00000000e+00],
        [2.71738575e-03, 2.39413914e+04, 6.30549922e+03, ...,
         1.46524884e+00, 6.71052823e-01, 4.89833757e+00],
        ...,
        [2.54381865e-03, 2.19837606e+04, 5.95692793e+03, ...,
         0.00000000e+00, 3.39133624e+01, 2.16510384e+01],
        [2.61253140e-03, 2.34605086e+04, 4.05115425e+03, ...,
         0.00000000e+00, 4.24551867e+01, 1.36400978e+01],
        [2.62450052e-03, 2.21588953e+04, 4.67899971e+03, ...,
         0.00000000e+00, 1.62856158e+01, 0.00000000e+00]]])
Coordinates:
  * time     (time) datetime64[ns] 2016-06-07T16:21:30
  * azimuth  (azimuth) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
  * range    (range) float64 29.98 89.94 149.9 ... 1.982e+04 1.988e+04 1.994e+04
Attributes: (12/17)
    units:                mm**6/m**3
    standard_name:        linear_equivalent_reflectivity_factor
    short_name:           Z
    Convections:          CF-1.0
    History:              2016-06-08 09:32 generated from orignal data
    Sample_Frequency_Hz:  5000000
    ...                   ...
    institution:          Meteorological Institute, Hamburg, Germany
    title:                LAWR data
    latitude:             53.56833
    longitude:            9.973997
    corr_offset:          29.5364
    tags:                 neg-lin-refl-repl-0.0;noise-filtered;beam-expansion...
Determined noiselevel: 3.714e-07

We notice that the logarithmic reflectivity was translated into a linear reflectivity with our noise remover. Thus, we need to retranslate our linear reflectivity to a logarithmic reflectivity.

reflectivity = reflectivity.lawr.to_dbz()

3. Filter out clutter

After we substracted a determined noise level from our reflectivity field, we want to remove clutter. Based on multiple application gradient-based filters andoptical filters, remove_clutter removes clutter from our reflectivity.

reflectivity, cluttermap = transform_funcs.remove_clutter_lawr(reflectivity)

print(reflectivity)
print(cluttermap)

<xarray.DataArray 'dbz' (time: 1, azimuth: 360, range: 333)>
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]])
Coordinates:
  * time     (time) datetime64[ns] 2016-06-07T16:21:30
  * azimuth  (azimuth) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
  * range    (range) float64 29.98 89.94 149.9 ... 1.982e+04 1.988e+04 1.994e+04
Attributes: (12/17)
    units:                dBZ
    standard_name:        equivalent_reflectivity_factor
    short_name:           dBZ
    Convections:          CF-1.0
    History:              2016-06-08 09:32 generated from orignal data
    Sample_Frequency_Hz:  5000000
    ...                   ...
    institution:          Meteorological Institute, Hamburg, Germany
    title:                LAWR data
    latitude:             53.56833
    longitude:            9.973997
    corr_offset:          29.5364
    tags:                 neg-lin-refl-repl-0.0;noise-filtered;beam-expansion...
ClutterMap(1*TDBZFilter, 1*SPINFilter, 1*SPKFilterOne, 1*SPKFilterTwo, 1*RINGFilterOne, 1*RINGFilterTwo, 1*SpeckleFilterW33, 1*SpeckleFilterW35, 1*SpeckleFilterW55, 1*SpeckleFilterW57, 1*SpeckleFilterW77)

4. Interpolate missing values

We can see that there are missing values within our reflectivity array, which need to be interpolated. For this, we will use interpolate_missing from our filter functions. This remaps existing values to missing values such that no missing value remains. If a criterion is exceeded, some missing values are filled with no rain to avoid interpolation artifacts.

reflectivity, remapper = transform_funcs.interpolate_missing(reflectivity)

print(reflectivity)
print(remapper)

/home/docs/checkouts/readthedocs.org/user_builds/pylawr/conda/stable/lib/python3.11/site-packages/xarray/core/coordinates.py:41: FutureWarning: Updating MultiIndexed coordinate 'grid_cell' would corrupt indices for other variables: ['azimuth', 'range']. This will raise an error in the future. Use `.drop_vars({'range', 'azimuth', 'grid_cell'})` before assigning new coordinate values.
  self.update({key: value})
/home/docs/checkouts/readthedocs.org/user_builds/pylawr/conda/stable/lib/python3.11/site-packages/xarray/core/coordinates.py:41: FutureWarning: Updating MultiIndexed coordinate 'grid_cell' would corrupt indices for other variables: ['grid_lat', 'grid_lon']. This will raise an error in the future. Use `.drop_vars({'grid_lat', 'grid_lon', 'grid_cell'})` before assigning new coordinate values.
  self.update({key: value})
<xarray.DataArray 'dbz' (time: 1, azimuth: 360, range: 333)>
array([[[-32.5, -32.5, -32.5, ..., -32.5, -32.5, -32.5],
        [-32.5, -32.5, -32.5, ..., -32.5, -32.5, -32.5],
        [-32.5, -32.5, -32.5, ..., -32.5, -32.5, -32.5],
        ...,
        [-32.5, -32.5, -32.5, ..., -32.5, -32.5, -32.5],
        [-32.5, -32.5, -32.5, ..., -32.5, -32.5, -32.5],
        [-32.5, -32.5, -32.5, ..., -32.5, -32.5, -32.5]]])
Coordinates:
  * azimuth  (azimuth) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
  * range    (range) float64 29.98 89.94 149.9 ... 1.982e+04 1.988e+04 1.994e+04
  * time     (time) datetime64[ns] 2016-06-07T16:21:30
Attributes: (12/17)
    units:                dBZ
    standard_name:        equivalent_reflectivity_factor
    short_name:           dBZ
    Convections:          CF-1.0
    History:              2016-06-08 09:32 generated from orignal data
    Sample_Frequency_Hz:  5000000
    ...                   ...
    institution:          Meteorological Institute, Hamburg, Germany
    title:                LAWR data
    latitude:             53.56833
    longitude:            9.973997
    corr_offset:          29.5364
    tags:                 neg-lin-refl-repl-0.0;noise-filtered;beam-expansion...
OrdinaryKriging(n_neighbors=10)

5. Plot the data

The missing values are interpolated with a nearest neighbor interpolation, where the five nearest neighbors are averaged. After we interpolated our data, we can plot the data.

To plot the data, we need to translate our reflectivity to rain rates. This translation is given by a Z-R relationship, where the scale and ploynomial degree factor can be set. All values without rain are set to NaN, because they will be plotted transparent.

rain_rate = reflectivity.lawr.to_z().lawr.zr_convert(a=256, b=1.42)
rain_rate = rain_rate.where(rain_rate >= 0.1)
rain_rate = rain_rate.lawr.set_grid_coordinates(grid)

After we translated the rain rate, we can use this array as input array to plot_rain_clutter function to plot this rain rate together with our cluttermap.

plotter = plot_funcs.plot_rain_clutter(rain_rate,
                                       cluttermap=cluttermap,
                                       plot_path='/tmp/quickstart_functional_api.png')
plotter.show()

This creates a default plot with an information header, a map with the rainfall rates and a corresponding colorbar.

After we plotted the rain rate and cluttermap, we finished this quickstart tutorial. If you want to get a deeper insight into this radar software, you can also check-out the quickstart tutorial for our low-level api.