.. note::
    :class: sphx-glr-download-link-note

    Click :ref:`here <sphx_glr_download_auto_examples_radar_qpe.py>` to download the full example code
.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_radar_qpe.py:


Radar QPE (Quantitative Precipitation Estimate)
===============================================
This example demonstrates how to generate radar QPE



.. code-block:: python


    import os
    import xarray as xr
    import matplotlib
    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    from pyresample import utils

    from matplotlib.colors import BoundaryNorm

    from swirlspy.rad.iris import read_iris
    from swirlspy.qpe.utils import dbz2rr, rr2rf, timestamps_ending, locate_file

    plt.switch_backend('agg')

    FILE_LOC = os.getcwd()
    THIS_DIR = FILE_LOC+'/../tests/int'
    os.chdir(THIS_DIR)





.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none

    ## You are using the Python ARM Radar Toolkit (Py-ART), an open source
    ## library for working with weather radar data. Py-ART is partly
    ## supported by the U.S. Department of Energy as part of the Atmospheric
    ## Radiation Measurement (ARM) Climate Research Facility, an Office of
    ## Science user facility.
    ##
    ## If you use this software to prepare a publication, please cite:
    ##
    ##     JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119


Read IRIS RAW data files into xarray
---------------------------------------------------

1. Generate timestamps with timestamps_ending()
2. Locate radar data files with locate_file()
3. Define target grid as a pyresample AreaDefinition
4. Read radar data files with read_iris()
5. Convert reflectivity in dBZ to rainrate in mm/h with dbz2rr()
6. Convert rainrate to rainfall in 6 mins with rr2rf()
7. Concatenate all xarrays into one with xarray.concat()
8. Compute the sum along the time dimension



.. code-block:: python


    dir = THIS_DIR+"/../samples/iris"
    located_files = []
    timestamps = timestamps_ending(pd.Timestamp(2018, 8, 22, 19, 00))
    for timestamp in timestamps:
        located_files.append(locate_file(dir, timestamp))

    area_id = "hk1980_250km"
    description = (
        "A 250 m resolution rectangular grid "
        "centred at HKO and extending to 250 km "
        "in each direction in HK1980 easting/northing coordinates"
        )
    proj_id = 'hk1980'
    projection = (
        '+proj=tmerc +lat_0=22.31213333333334 '
        '+lon_0=114.1785555555556 +k=1 +x_0=836694.05 '
        '+y_0=819069.8 +ellps=intl +towgs84=-162.619,-276.959,'
        '-161.764,0.067753,-2.24365,-1.15883,-1.09425 +units=m '
        '+no_defs'
        )
    x_size = 500
    y_size = 500
    area_extent = (587000, 569000, 1087000, 1069000)
    area_def = utils.get_area_def(
        area_id, description, proj_id, projection, x_size, y_size, area_extent
    )

    rainfalls = []
    for filename in located_files:
        reflectivity = read_iris(
            filename,
            area_def=area_def,
            coord_label=['easting', 'northing'],
            radar_radius=250000
        )
        rainrate = dbz2rr(reflectivity, a=58.53, b=1.56)
        rainfall = rr2rf(rainrate, time_shift=pd.Timedelta(minutes=6))
        rainfalls.append(rainfall)

    rainfall_cat = xr.concat(rainfalls, dim='time')
    rainfall_hourly = rainfall_cat.sum(dim='time')







Generate rainfall map
---------------------------------------------------
Define colour scale and format of the rainfall map plot
then generate it with xarray.plot()



.. code-block:: python


    levels = [
        0, 0.5, 2, 5, 10, 20,
        30, 40, 50, 70, 100, 150,
        200, 300, 400, 500, 600, 700
    ]

    cmap = matplotlib.colors.ListedColormap([
        '#ffffff', '#9bf7f7', '#00ffff', '#00d5cc', '#00bd3d', '#2fd646',
        '#9de843', '#ffdd41', '#ffac33', '#ff621e', '#d23211', '#9d0063',
        '#e300ae', '#ff00ce', '#ff57da', '#ff8de6', '#ffe4fd'
    ])

    norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
    plt.figure(figsize=(15, 10))
    plt.axis('equal')
    rainfall_hourly.plot(cmap=cmap, norm=norm)
    plt.savefig("radar_qpe.png")



.. image:: /auto_examples/images/sphx_glr_radar_qpe_001.png
    :class: sphx-glr-single-img




**Total running time of the script:** ( 0 minutes  19.456 seconds)


.. _sphx_glr_download_auto_examples_radar_qpe.py:


.. only :: html

 .. container:: sphx-glr-footer
    :class: sphx-glr-footer-example



  .. container:: sphx-glr-download

     :download:`Download Python source code: radar_qpe.py <radar_qpe.py>`



  .. container:: sphx-glr-download

     :download:`Download Jupyter notebook: radar_qpe.ipynb <radar_qpe.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.readthedocs.io>`_