Interpolate data on cartesian coordinates to polar coordinates

import numpy as np
import matplotlib.pyplot as plt
import wradlib as wrl
import wradlib_data
import warnings
import xarray as xr
import cmweather

warnings.filterwarnings("ignore")
try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    plt.ion()

/home/docs/checkouts/readthedocs.org/user_builds/wradlib/conda/latest/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Read a cartesian data set

filename = wradlib_data.DATASETS.fetch("geo/bonn_new.tif")
print(filename)

Downloading file 'geo/bonn_new.tif' from 'https://github.com/wradlib/wradlib-data/raw/main/data/geo/bonn_new.tif' to '/home/docs/.cache/wradlib-data'.

/home/docs/.cache/wradlib-data/geo/bonn_new.tif

dem_bonn = xr.open_dataset(filename, engine="rasterio")

Inspect the data set a little

display(dem_bonn)

<xarray.Dataset> Size: 69MB
Dimensions:      (band: 1, x: 4800, y: 3600)
Coordinates:
  * band         (band) int64 8B 1
  * x            (x) float64 38kB 5.0 5.001 5.002 5.003 ... 8.998 8.999 9.0
  * y            (y) float64 29kB 52.0 52.0 52.0 52.0 ... 49.0 49.0 49.0 49.0
    spatial_ref  int64 8B ...
Data variables:
    band_data    (band, y, x) float32 69MB ...

xarray.Dataset

• Dimensions:
  • band: 1
  • x: 4800
  • y: 3600
• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • x
    (x)
    float64
    5.0 5.001 5.002 ... 8.998 8.999 9.0

    array([5.000417, 5.00125 , 5.002083, ..., 8.997917, 8.99875 , 8.999583],
          shape=(4800,))

  • y
    (y)
    float64
    52.0 52.0 52.0 ... 49.0 49.0 49.0

    array([51.999583, 51.99875 , 51.997917, ..., 49.002083, 49.00125 , 49.000417],
          shape=(3600,))

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    5.0 0.0008333333333333334 0.0 52.0 0.0 -0.0008333333333333334

    [1 values with dtype=int64]

• Data variables: (1)
  • band_data
    (band, y, x)
    float32
    ...

    AREA_OR_POINT :

    Area

    [17280000 values with dtype=float32]

Extract dem_bonn band and coarsen grid. We do this here to prevent memory issues when running this on CI. When applying this to your workflows, just use your normal grid resolution.

 band = dem_bonn.coarsen(x=10, y=10, boundary="trim").mean()["band_data"].isel(band=0)
 display(band)

<xarray.DataArray 'band_data' (y: 360, x: 480)> Size: 691kB
array([[-9.0000e-01, -1.6100e+00, -1.5000e+00, ...,  1.5455e+02,
         1.6103e+02,  1.7505e+02],
       [-4.5000e-01, -3.6000e-01,  1.8000e-01, ...,  1.3796e+02,
         1.4295e+02,  1.4829e+02],
       [ 1.0300e+00,  8.7000e-01,  4.4000e-01, ...,  1.5510e+02,
         1.4532e+02,  1.6954e+02],
       ...,
       [ 1.6030e+02,  1.6413e+02,  1.6526e+02, ...,  2.8768e+02,
         2.9277e+02,  3.0780e+02],
       [ 1.5684e+02,  1.5749e+02,  1.6067e+02, ...,  3.0840e+02,
         2.9680e+02,  2.7706e+02],
       [ 1.5536e+02,  1.5613e+02,  1.6143e+02, ...,  3.8239e+02,
         4.3090e+02,  3.6198e+02]], shape=(360, 480), dtype=float32)
Coordinates:
  * y            (y) float64 3kB 52.0 51.99 51.98 51.97 ... 49.02 49.01 49.0
  * x            (x) float64 4kB 5.004 5.013 5.021 5.029 ... 8.979 8.988 8.996
    band         int64 8B 1
    spatial_ref  int64 8B ...
Attributes:
    AREA_OR_POINT:  Area

xarray.DataArray

'band_data'

• y: 360
• x: 480

• -0.9 -1.61 -1.5 -1.13 -0.88 -0.4 ... 334.3 357.5 382.4 430.9 362.0

  array([[-9.0000e-01, -1.6100e+00, -1.5000e+00, ...,  1.5455e+02,
           1.6103e+02,  1.7505e+02],
         [-4.5000e-01, -3.6000e-01,  1.8000e-01, ...,  1.3796e+02,
           1.4295e+02,  1.4829e+02],
         [ 1.0300e+00,  8.7000e-01,  4.4000e-01, ...,  1.5510e+02,
           1.4532e+02,  1.6954e+02],
         ...,
         [ 1.6030e+02,  1.6413e+02,  1.6526e+02, ...,  2.8768e+02,
           2.9277e+02,  3.0780e+02],
         [ 1.5684e+02,  1.5749e+02,  1.6067e+02, ...,  3.0840e+02,
           2.9680e+02,  2.7706e+02],
         [ 1.5536e+02,  1.5613e+02,  1.6143e+02, ...,  3.8239e+02,
           4.3090e+02,  3.6198e+02]], shape=(360, 480), dtype=float32)

• Coordinates: (4)
  • y
    (y)
    float64
    52.0 51.99 51.98 ... 49.01 49.0

    array([51.995833, 51.9875  , 51.979167, ..., 49.020833, 49.0125  , 49.004167],
          shape=(360,))

  • x
    (x)
    float64
    5.004 5.013 5.021 ... 8.988 8.996

    array([5.004167, 5.0125  , 5.020833, ..., 8.979167, 8.9875  , 8.995833],
          shape=(480,))

  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    5.0 0.0008333333333333334 0.0 52.0 0.0 -0.0008333333333333334

    [1 values with dtype=int64]

• Attributes: (1)

  AREA_OR_POINT :

  Area

Read a polar data set

In order to be on the same coordinates, we georeference our radar sweep accordingly.

fname = wradlib_data.DATASETS.fetch("hdf5/2014-08-10--182000.ppi.mvol")
swp = xr.open_dataset(fname, engine="gamic", group="sweep_0")
swp = swp.wrl.georef.georeference(crs=dem_bonn.spatial_ref)
swp = swp.set_coords("sweep_mode")
display(swp)

<xarray.Dataset> Size: 33MB
Dimensions:            (azimuth: 360, range: 1000)
Coordinates: (12/15)
  * azimuth            (azimuth) float64 3kB 0.5054 1.522 2.527 ... 358.5 359.5
    elevation          (azimuth) float64 3kB 1.505 1.505 1.505 ... 1.505 1.505
    time               (azimuth) datetime64[ns] 3kB ...
  * range              (range) float32 4kB 50.0 150.0 ... 9.985e+04 9.995e+04
    x                  (azimuth, range) float64 3MB 7.072 7.072 ... 7.059 7.059
    y                  (azimuth, range) float64 3MB 50.73 50.73 ... 51.63 51.63
    ...                 ...
    bins               (azimuth, range) float32 1MB 50.0 150.0 ... 9.995e+04
    sweep_mode         <U20 80B 'azimuth_surveillance'
    longitude          float64 8B 7.072
    latitude           float64 8B 50.73
    altitude           float64 8B 99.5
    crs_wkt            int64 8B 0
Data variables: (12/17)
    KDP                (azimuth, range) float32 1MB ...
    PHIDP              (azimuth, range) float32 1MB ...
    DBZH               (azimuth, range) float32 1MB ...
    DBZV               (azimuth, range) float32 1MB ...
    RHOHV              (azimuth, range) float32 1MB ...
    DBTH               (azimuth, range) float32 1MB ...
    ...                 ...
    ZDR                (azimuth, range) float32 1MB ...
    sweep_number       int64 8B ...
    prt_mode           <U7 28B ...
    follow_mode        <U7 28B ...
    sweep_fixed_angle  float64 8B ...
    nyquist_velocity   object 8B ...
Attributes:
    source:       gamic
    pulse_width:  2

xarray.Dataset

• Dimensions:
  • azimuth: 360
  • range: 1000
• Coordinates: (15)
  • azimuth
    (azimuth)
    float64
    0.5054 1.522 2.527 ... 358.5 359.5

    standard_name :

    ray_azimuth_angle

    long_name :

    azimuth_angle_from_true_north

    units :

    degrees

    axis :

    radial_azimuth_coordinate

    array([  0.505371,   1.521606,   2.526855, ..., 357.50885 , 358.514099,
           359.502869], shape=(360,))

  • elevation
    (azimuth)
    float64
    1.505 1.505 1.505 ... 1.505 1.505

    standard_name :

    ray_elevation_angle

    long_name :

    elevation_angle_from_horizontal_plane

    units :

    degrees

    axis :

    radial_elevation_coordinate

    array([1.505127, 1.505127, 1.505127, ..., 1.505127, 1.505127, 1.505127],
          shape=(360,))

  • time
    (azimuth)
    datetime64[ns]
    ...

    standard_name :

    time

    [360 values with dtype=datetime64[ns]]

  • range
    (range)
    float32
    50.0 150.0 ... 9.985e+04 9.995e+04

    units :

    meters

    standard_name :

    projection_range_coordinate

    long_name :

    range_to_measurement_volume

    axis :

    radial_range_coordinate

    meters_between_gates :

    100.0

    spacing_is_constant :

    true

    meters_to_center_of_first_gate :

    50.0

    array([5.000e+01, 1.500e+02, 2.500e+02, ..., 9.975e+04, 9.985e+04, 9.995e+04],
          shape=(1000,), dtype=float32)

  • x
    (azimuth, range)
    float64
    7.072 7.072 7.072 ... 7.059 7.059

    standard_name :

    longitude

    long_name :

    longitude coordinate

    units :

    degrees_east

    axis :

    X

    array([[7.07166924, 7.07168173, 7.07169422, ..., 7.0843588 , 7.08437177,
            7.08438475],
           [7.0716818 , 7.0717194 , 7.071757  , ..., 7.10988417, 7.10992322,
            7.10996228],
           [7.07169421, 7.07175664, 7.07181906, ..., 7.13512108, 7.13518592,
            7.13525076],
           ...,
           [7.07163223, 7.07157069, 7.07150914, ..., 7.00910101, 7.00903708,
            7.00897316],
           [7.07164464, 7.07160793, 7.07157121, ..., 7.0343385 , 7.03430036,
            7.03426223],
           [7.07165686, 7.07164457, 7.07163229, ..., 7.05917419, 7.05916143,
            7.05914867]], shape=(360, 1000))

  • y
    (azimuth, range)
    float64
    50.73 50.73 50.73 ... 51.63 51.63

    standard_name :

    latitude

    long_name :

    latitude coordinate

    units :

    degrees_north

    axis :

    Y

    array([[50.73096929, 50.73186786, 50.73276644, ..., 51.62646241,
            51.62736017, 51.62825793],
           [50.73096915, 50.73186744, 50.73276573, ..., 51.62617586,
            51.62707332, 51.62797079],
           [50.73096887, 50.73186661, 50.73276434, ..., 51.62560981,
            51.6265067 , 51.62740359],
           ...,
           [50.73096888, 50.73186664, 50.7327644 , ..., 51.62563473,
            51.62653164, 51.62742856],
           [50.73096915, 50.73186746, 50.73276577, ..., 51.6261908 ,
            51.62708828, 51.62798576],
           [50.73096929, 50.73186787, 50.73276644, ..., 51.62646356,
            51.62736132, 51.62825908]], shape=(360, 1000))

  • z
    (azimuth, range)
    float64
    100.8 103.4 ... 3.309e+03 3.313e+03

    standard_name :

    height_above_ground

    units :

    degrees_east

    array([[ 100.81346717,  103.44128457,  106.07027939, ..., 3305.14502024,
            3308.94609379, 3312.74837169],
           [ 100.81346717,  103.44128457,  106.07027939, ..., 3305.14502024,
            3308.94609379, 3312.74837169],
           [ 100.81346717,  103.44128457,  106.07027939, ..., 3305.14502024,
            3308.94609379, 3312.74837169],
           ...,
           [ 100.81346717,  103.44128457,  106.07027939, ..., 3305.14502024,
            3308.94609379, 3312.74837169],
           [ 100.81346717,  103.44128457,  106.07027939, ..., 3305.14502024,
            3308.94609379, 3312.74837169],
           [ 100.81346717,  103.44128457,  106.07027939, ..., 3305.14502024,
            3308.94609379, 3312.74837169]], shape=(360, 1000))

  • gr
    (azimuth, range)
    float64
    49.98 149.9 ... 9.978e+04 9.988e+04

    standard_name :

    distance_from_radar

    units :

    degrees_east

    array([[4.99828990e+01, 1.49947163e+02, 2.49911365e+02, ...,
            9.96790594e+04, 9.97789482e+04, 9.98788369e+04],
           [4.99829007e+01, 1.49947165e+02, 2.49911367e+02, ...,
            9.96790594e+04, 9.97789482e+04, 9.98788369e+04],
           [4.99829022e+01, 1.49947166e+02, 2.49911368e+02, ...,
            9.96790594e+04, 9.97789482e+04, 9.98788369e+04],
           ...,
           [4.99828926e+01, 1.49947157e+02, 2.49911359e+02, ...,
            9.96790594e+04, 9.97789482e+04, 9.98788369e+04],
           [4.99828950e+01, 1.49947159e+02, 2.49911361e+02, ...,
            9.96790594e+04, 9.97789482e+04, 9.98788369e+04],
           [4.99828971e+01, 1.49947161e+02, 2.49911363e+02, ...,
            9.96790594e+04, 9.97789482e+04, 9.98788369e+04]],
          shape=(360, 1000))

  • rays
    (azimuth, range)
    float64
    0.5054 0.5054 ... 359.5 359.5

    standard_name :

    ray_azimuth_angle

    long_name :

    azimuth_angle_from_true_north

    units :

    degrees

    axis :

    radial_azimuth_coordinate

    array([[  0.50537109,   0.50537109,   0.50537109, ...,   0.50537109,
              0.50537109,   0.50537109],
           [  1.52160645,   1.52160645,   1.52160645, ...,   1.52160645,
              1.52160645,   1.52160645],
           [  2.52685547,   2.52685547,   2.52685547, ...,   2.52685547,
              2.52685547,   2.52685547],
           ...,
           [357.5088501 , 357.5088501 , 357.5088501 , ..., 357.5088501 ,
            357.5088501 , 357.5088501 ],
           [358.51409912, 358.51409912, 358.51409912, ..., 358.51409912,
            358.51409912, 358.51409912],
           [359.50286865, 359.50286865, 359.50286865, ..., 359.50286865,
            359.50286865, 359.50286865]], shape=(360, 1000))

  • bins
    (azimuth, range)
    float32
    50.0 150.0 ... 9.985e+04 9.995e+04

    units :

    meters

    standard_name :

    projection_range_coordinate

    long_name :

    range_to_measurement_volume

    axis :

    radial_range_coordinate

    meters_between_gates :

    100.0

    spacing_is_constant :

    true

    meters_to_center_of_first_gate :

    50.0

    array([[5.000e+01, 1.500e+02, 2.500e+02, ..., 9.975e+04, 9.985e+04,
            9.995e+04],
           [5.000e+01, 1.500e+02, 2.500e+02, ..., 9.975e+04, 9.985e+04,
            9.995e+04],
           [5.000e+01, 1.500e+02, 2.500e+02, ..., 9.975e+04, 9.985e+04,
            9.995e+04],
           ...,
           [5.000e+01, 1.500e+02, 2.500e+02, ..., 9.975e+04, 9.985e+04,
            9.995e+04],
           [5.000e+01, 1.500e+02, 2.500e+02, ..., 9.975e+04, 9.985e+04,
            9.995e+04],
           [5.000e+01, 1.500e+02, 2.500e+02, ..., 9.975e+04, 9.985e+04,
            9.995e+04]], shape=(360, 1000), dtype=float32)

  • sweep_mode
    ()
    <U20
    'azimuth_surveillance'

    array('azimuth_surveillance', dtype='<U20')

  • longitude
    ()
    float64
    7.072

    long_name :

    longitude

    units :

    degrees_east

    standard_name :

    longitude

    array(7.071663)

  • latitude
    ()
    float64
    50.73

    long_name :

    latitude

    units :

    degrees_north

    positive :

    up

    standard_name :

    latitude

    array(50.73052)

  • altitude
    ()
    float64
    99.5

    long_name :

    altitude

    units :

    meters

    standard_name :

    altitude

    array(99.5)

  • crs_wkt
    ()
    int64
    0

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    array(0)

• Data variables: (17)
  • KDP
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Specific differential phase HV

    units :

    degrees per kilometer

    standard_name :

    radar_specific_differential_phase_hv

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • PHIDP
    (azimuth, range)
    float32
    ...

    format :

    UV16

    long_name :

    Differential phase HV

    units :

    degrees

    standard_name :

    radar_differential_phase_hv

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • DBZH
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Equivalent reflectivity factor H

    units :

    dBZ

    standard_name :

    radar_equivalent_reflectivity_factor_h

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • DBZV
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Equivalent reflectivity factor V

    units :

    dBZ

    standard_name :

    radar_equivalent_reflectivity_factor_v

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • RHOHV
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Correlation coefficient HV

    units :

    unitless

    standard_name :

    radar_correlation_coefficient_hv

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • DBTH
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Total power H (uncorrected reflectivity)

    units :

    dBZ

    standard_name :

    radar_equivalent_reflectivity_factor_h

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • DBTV
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Total power V (uncorrected reflectivity)

    units :

    dBZ

    standard_name :

    radar_equivalent_reflectivity_factor_v

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • VRADH
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Radial velocity of scatterers away from instrument H

    units :

    meters per seconds

    standard_name :

    radial_velocity_of_scatterers_away_from_instrument_h

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • VRADV
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Radial velocity of scatterers away from instrument V

    units :

    meters per second

    standard_name :

    radial_velocity_of_scatterers_away_from_instrument_v

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • WRADH
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Doppler spectrum width H

    units :

    meters per seconds

    standard_name :

    radar_doppler_spectrum_width_h

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • WRADV
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Doppler spectrum width V

    units :

    meters per second

    standard_name :

    radar_doppler_spectrum_width_v

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • ZDR
    (azimuth, range)
    float32
    ...

    format :

    UV8

    long_name :

    Log differential reflectivity H/V

    units :

    dB

    standard_name :

    radar_differential_reflectivity_hv

    _Undetect :

    0.0

    [360000 values with dtype=float32]

  • sweep_number
    ()
    int64
    ...

    [1 values with dtype=int64]

  • prt_mode
    ()
    <U7
    ...

    [1 values with dtype=<U7]

  • follow_mode
    ()
    <U7
    ...

    [1 values with dtype=<U7]

  • sweep_fixed_angle
    ()
    float64
    ...

    [1 values with dtype=float64]

  • nyquist_velocity
    ()
    object
    ...

    standard_name :

    nyquist_velocity

    units :

    m s-1

    [1 values with dtype=object]

• Attributes: (2)

  source :

  gamic

  pulse_width :

  2

griddata

For details see \(\omega radlib\)’s central interpolation API entrypoint wradlib.ipol.IpolMethods.interpolate. This is slow for large input arrays. Please check with scipy.interpolate.griddata for kwarg distribution.

band_xy = band.chunk()
swp = swp.isel(range=slice(0, 100))
band_polar_nearest = band_xy.wrl.ipol.interpolate(swp, method="griddata_nearest")
band_polar_linear = band_xy.wrl.ipol.interpolate(swp, method="griddata_linear")
band_polar_cubic = band_xy.wrl.ipol.interpolate(swp, method="griddata_cubic")

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 10))
cmap = "terrain"
vmin, vmax = 0.0, 250.0

band_crop = band_xy.wrl.util.crop(swp)
pm = band_crop.plot(ax=ax1, cmap=cmap, vmin=vmin, vmax=vmax)
ax1.set_aspect(band_crop.wrl.util.aspect())
ax1.set_title("Original DEM")

pm = band_polar_nearest.wrl.vis.plot(ax=ax2, cmap=cmap, vmin=vmin, vmax=vmax)
ax2.set_title("Nearest")
pm = band_polar_linear.wrl.vis.plot(ax=ax3, cmap=cmap, vmin=vmin, vmax=vmax)
ax3.set_title("Linear")
pm = band_polar_cubic.wrl.vis.plot(ax=ax4, cmap=cmap, vmin=vmin, vmax=vmax)
ax4.set_title("Cubic")

plt.tight_layout()

map_coordinates

For details see \(\omega radlib\)’s central interpolation API entrypoint wradlib.ipol.IpolMethods.interpolate. Please check with scipy.ndimage.map_coordinates for kwarg distribution.

band_xy = band.chunk()
swp = swp.isel(range=slice(0,100))
band_polar_nearest = band_xy.wrl.ipol.interpolate(swp, method="map_coordinates", order=0)
band_polar_linear = band_xy.wrl.ipol.interpolate(swp, method="map_coordinates", order=1)
band_polar_quadratic = band_xy.wrl.ipol.interpolate(swp, method="map_coordinates", order=2)
band_polar_cubic = band_xy.wrl.ipol.interpolate(swp, method="map_coordinates", order=3)
band_polar_quartic = band_xy.wrl.ipol.interpolate(swp, method="map_coordinates", order=4)
band_polar_quintic = band_xy.wrl.ipol.interpolate(swp, method="map_coordinates", order=5)

fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(12, 15))

pm = band_polar_nearest.wrl.vis.plot(ax=ax1, cmap=cmap, vmin=vmin, vmax=vmax)
ax1.set_title("Nearest")
pm = band_polar_linear.wrl.vis.plot(ax=ax2, cmap=cmap, vmin=vmin, vmax=vmax)
ax2.set_title("Linear")
pm = band_polar_quadratic.wrl.vis.plot(ax=ax3, cmap=cmap, vmin=vmin, vmax=vmax)
ax3.set_title("Quadratic")
pm = band_polar_cubic.wrl.vis.plot(ax=ax4, cmap=cmap, vmin=vmin, vmax=vmax)
ax4.set_title("Cubic")
pm = band_polar_quartic.wrl.vis.plot(ax=ax5, cmap=cmap, vmin=vmin, vmax=vmax)
ax5.set_title("Quartic")
pm = band_polar_quintic.wrl.vis.plot(ax=ax6, cmap=cmap, vmin=vmin, vmax=vmax)
ax6.set_title("Quintic")

plt.tight_layout()