.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/qpe_manila.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_qpe_manila.py: QPE (Manila) ==================================================== This example demonstrates how to perform QPE, using raingauge data from Manila and radar data from Subic. .. GENERATED FROM PYTHON SOURCE LINES 10-13 Definitions -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 13-44 .. code-block:: default import os import time import numpy as np import pandas as pd import copy import xarray import scipy import pyproj import matplotlib from matplotlib.colors import BoundaryNorm import matplotlib.pyplot as plt from PIL import Image from pyresample import utils from cartopy.io import shapereader import cartopy.feature as cfeature import cartopy.crs as ccrs import cartopy.io.shapereader as shpreader from swirlspy.obs import Rain from swirlspy.rad.uf_ph import read_uf_ph from swirlspy.core.resample import grid_resample from swirlspy.qpe.rfmap import rg_interpolate, comp_qpe, show_raingauge from swirlspy.qpe.utils import timestamps_ending, locate_file, dbz2rr, rr2rf, \ temporal_interp, acc plt.switch_backend('agg') THIS_DIR = os.getcwd() os.chdir(THIS_DIR) .. GENERATED FROM PYTHON SOURCE LINES 45-52 Initialising ----------------------------------------------------------------- This section demonstrates extracting raingauge and radar data. Step 1: Defining an end-time for accumulating rainfall. .. GENERATED FROM PYTHON SOURCE LINES 52-56 .. code-block:: default acctime = pd.Timestamp('20180811112000').floor('min') acctime_str = acctime.strftime('%Y%m%d%H%M') .. GENERATED FROM PYTHON SOURCE LINES 57-59 Step 2: Setting up directories for raingauge and radar files. .. GENERATED FROM PYTHON SOURCE LINES 59-63 .. code-block:: default rad_dir = THIS_DIR + '/../tests/samples/uf_ph/sub/' rg_dir = THIS_DIR + '/../tests/samples/rfmap/' .. GENERATED FROM PYTHON SOURCE LINES 64-66 Step 3: Generating timestamps and pattern for both radar and raingauge files. .. GENERATED FROM PYTHON SOURCE LINES 66-92 .. code-block:: default # Timestamps of raingauges rg_timestrings = timestamps_ending( acctime + pd.Timedelta(minutes=10), duration=pd.Timedelta(hours=1), interval=pd.Timedelta(minutes=10) ) # Raingauge pattern rg_pattern = ['rf60m_20'+ts for ts in rg_timestrings] # Finding time nearest radar file # to accumulation end time minute = acctime.minute nearest_6min = acctime.minute // 10 * 10 nearest_rad_timestamp = pd.Timestamp( acctime_str[:-2]+f'{nearest_6min:02}' ) rad_timestrings = timestamps_ending( nearest_rad_timestamp, duration=pd.Timedelta(hours=1), interval=pd.Timedelta(minutes=10) ) .. GENERATED FROM PYTHON SOURCE LINES 93-95 Step 4: Extracting raingauge and radar files from their respective directories. .. GENERATED FROM PYTHON SOURCE LINES 95-106 .. code-block:: default located_rg_files = [] for pat in rg_pattern: located_rg_files.append(locate_file(rg_dir, pat)) located_radar_files = [] for ts in rad_timestrings: located_radar_files.append(locate_file(rad_dir, ts)) .. GENERATED FROM PYTHON SOURCE LINES 107-109 Step 5: Read data from raingauge files into a Rain object. Coordinates are geodetic, following that in the files. .. GENERATED FROM PYTHON SOURCE LINES 109-117 .. code-block:: default rg_object_geodetic = Rain( located_rg_files, 'WGS84', duration=pd.Timedelta(minutes=5), NAN=[3276.7, 32767] ) .. GENERATED FROM PYTHON SOURCE LINES 118-119 Step 6: Define the target grid as a pyresample AreaDefinition. .. GENERATED FROM PYTHON SOURCE LINES 119-138 .. code-block:: default # Defining target grid area_id = "epsg3123_240km" description = ("A 240 m resolution rectangular grid " "centred at Subic RADAR and extending to 240 km " "in each direction") proj_id = 'epsg3123' projection = ('+proj=tmerc +lat_0=0 ' '+lon_0=121 +k=0.99995 +x_0=500000 ' '+y_0=0 +ellps=clrk66 +towgs84=-127.62,-67.24,' '-47.04,-3.068,4.903,1.578,-1.06 +units=m ' '+no_defs') x_size = 1000 y_size = 1000 area_extent = (191376.04113, 1399386.68659, 671376.04113, 1879386.68659) area_def = utils.get_area_def( area_id, description, proj_id, projection, x_size, y_size, area_extent ) .. rst-class:: sphx-glr-script-out .. code-block:: none /tmp/build/docs/swirlspy/swirlspy/examples/qpe_manila.py:135: UserWarning: 'get_area_def' has moved, import it with 'from pyresample import get_area_def' area_id, description, proj_id, projection, x_size, y_size, area_extent /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() .. GENERATED FROM PYTHON SOURCE LINES 139-143 Step 7: Convert coordinates of raingauge object to desired projection. In this example, the desired projection is PRS92. This can be achieved by the .reproject() method of the Rain object. .. GENERATED FROM PYTHON SOURCE LINES 143-149 .. code-block:: default inProj = pyproj.Proj(init="epsg:4326") outProj = pyproj.Proj(area_def.proj_str) rg_object = rg_object_geodetic.reproject(inProj, outProj, "PRS92") .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs)) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:294: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 projstring = _prepare_from_string(" ".join((projstring, projkwargs))) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() .. GENERATED FROM PYTHON SOURCE LINES 150-151 Step 8: Read radar files into xarray.DataArrays using read_uf_ph(). .. GENERATED FROM PYTHON SOURCE LINES 151-164 .. code-block:: default reflec_list = [] for file in located_radar_files: reflec = read_uf_ph( file, area_def=area_def, coord_label=['easting', 'northing'], indicator='deg_ppi', elevation=0.5 ) reflec_list.append(reflec) reflectivity = xarray.concat(reflec_list, dim='time') .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs)) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:294: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 projstring = _prepare_from_string(" ".join((projstring, projkwargs))) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs)) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:294: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 projstring = _prepare_from_string(" ".join((projstring, projkwargs))) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs)) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:294: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 projstring = _prepare_from_string(" ".join((projstring, projkwargs))) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs)) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:294: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 projstring = _prepare_from_string(" ".join((projstring, projkwargs))) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs)) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:294: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 projstring = _prepare_from_string(" ".join((projstring, projkwargs))) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs)) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:294: FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 projstring = _prepare_from_string(" ".join((projstring, projkwargs))) /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() .. GENERATED FROM PYTHON SOURCE LINES 165-173 Accumulating and interpolating rainfall ----------------------------------------------------------------- Interpolate rainfall recorded by raingauges into the user-defined grid and accumulate radar rainfall over an hour after making the necessary adjustments. .. GENERATED FROM PYTHON SOURCE LINES 175-179 Step 1: Interpolate Rain object to user-defined grid. In this example, a multiquadric Radial Basis Function is used. .. GENERATED FROM PYTHON SOURCE LINES 179-186 .. code-block:: default interpolated_rg = rg_interpolate( rg_object, area_def, 'rbf', coord_label=['easting', 'northing'] ) .. GENERATED FROM PYTHON SOURCE LINES 187-190 Step 2: Convert to radar reflectivity to rainrates, convert rainrates to times of raingauges, and accumulate rainfalls every 10 minutes. .. GENERATED FROM PYTHON SOURCE LINES 190-203 .. code-block:: default rainrates = dbz2rr(reflectivity, a=300, b=1.4) # Convert time coordinates of rainrates from start time # to end time rainrates_time_endtime = rainrates.copy() rainrates_time_endtime.coords['time'] = \ [ pd.Timestamp(t) + pd.Timedelta(minutes=10) for t in rainrates.coords['time'].values ] rainfalls = rr2rf(rainrates_time_endtime, scan_duration=10) .. GENERATED FROM PYTHON SOURCE LINES 204-205 Step 3: Accumulate rainfall over an hour. .. GENERATED FROM PYTHON SOURCE LINES 205-212 .. code-block:: default acc_rf = acc( rainfalls, rg_object.end_time, acc_period=pd.Timedelta(minutes=60) ) .. GENERATED FROM PYTHON SOURCE LINES 213-218 Compositing rainfall ----------------------------------------------------------------- Perform compositing on radar and raingauge derived rainfall to obtain a composite QPE. .. GENERATED FROM PYTHON SOURCE LINES 218-226 .. code-block:: default comprf = comp_qpe( area_def, rg_object=rg_object, rg_interp=interpolated_rg, rad_rf=acc_rf ) .. GENERATED FROM PYTHON SOURCE LINES 227-232 Plotting --------------------------------------------------------------- Plot composited radar and raingauge rainfall. .. GENERATED FROM PYTHON SOURCE LINES 232-352 .. code-block:: default # Defining coastlines map_shape_file = os.path.join(THIS_DIR, "./../tests/samples/shape/rsmc") ocean_color = np.array([[[178, 208, 254]]], dtype=np.uint8) land_color = cfeature.COLORS['land'] coastline = cfeature.ShapelyFeature( list(shpreader.Reader(map_shape_file).geometries()), ccrs.PlateCarree() ) # Plotting function for neatness. def plot_map( da, rg_object, acctime, area_def, based='raingauge and radar', savepath='', area_extent=None ): """ A custom function for plotting a map. Parameters -------------- da: xarray.DataArray Contains data to be plotted. rg_object: Rain Contains raingauge data. acctime: pd.Timestamp Contains the endtime of the accumulation period. area_def: pyresample.geometry.AreaDefinition AreaDefinition of the grid. based: str Type of data plotted in the map. savepath: str Path to save the image to. area_extent: tuple Area extent (x0, x1, y0, y1) to be plotted. Defaults to None. """ # Defining the colour scheme 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) # Plotting axes plt.figure(figsize=(24, 21)) crs = area_def.to_cartopy_crs() ax = plt.axes(projection=crs) if area_extent is not None: ax.set_extent(area_extent, crs=crs) # Plot data quadmesh = da.plot( cmap=cmap, norm=norm, extend='max', cbar_kwargs={'ticks': levels, 'format': '%.3g'} ) # Adjusting size of colorbar cb = quadmesh.colorbar cb.ax.set_ylabel( da.attrs['long_name']+'['+da.attrs['units']+']', fontsize=28 ) cb.ax.tick_params(labelsize=24) # Setting labels ax.xaxis.set_visible(True) ax.yaxis.set_visible(True) for tick in ax.xaxis.get_major_ticks(): tick.label.set_fontsize(24) for tick in ax.yaxis.get_major_ticks(): tick.label.set_fontsize(24) ax.xaxis.label.set_size(28) ax.yaxis.label.set_size(28) # ocean ax.imshow(np.tile(ocean_color, [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180], zorder=-1) # coastline, color ax.add_feature(coastline, facecolor=land_color, edgecolor='none', zorder=0) # overlay coastline without color ax.add_feature(coastline, facecolor='none', edgecolor='gray', linewidth=0.5, zorder=3) # Show raingauge show_raingauge( rg_object, ax, show_value=True, color='red', markersize=20, fontsize=20 ) # Show title plt.title( (f"Last Hour Rainfall at {acctime.strftime('%H:%MH %d-%b-%Y')}\n" f"(based on {based} data)"), fontsize=32 ) plt.savefig(savepath) .. GENERATED FROM PYTHON SOURCE LINES 353-390 .. code-block:: default # Plotting maps r = 64000 proj_site = reflectivity.proj_site area_extent = ( proj_site[0]-r, proj_site[0]+r, proj_site[1]-r, proj_site[1]+r ) # Raingauge only plot_map( interpolated_rg, rg_object, acctime, area_def, based='raingauge', savepath=THIS_DIR+f'/../tests/outputs/raingauge_{acctime_str}.png', area_extent=area_extent ) # Radar only plot_map( acc_rf, rg_object, acctime, area_def, based='radar', savepath=THIS_DIR+f'/../tests/outputs/radar_{acctime_str}.png', area_extent=area_extent ) # Composite raingauge and radar plot_map( comprf, rg_object, acctime, area_def, savepath=THIS_DIR+f'/../tests/outputs/comp_{acctime_str}.png', area_extent=area_extent ) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_qpe_manila_001.png :alt: Last Hour Rainfall at 11:20H 11-Aug-2018 (based on raingauge data) :srcset: /auto_examples/images/sphx_glr_qpe_manila_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_qpe_manila_002.png :alt: Last Hour Rainfall at 11:20H 11-Aug-2018 (based on radar data) :srcset: /auto_examples/images/sphx_glr_qpe_manila_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_qpe_manila_003.png :alt: Last Hour Rainfall at 11:20H 11-Aug-2018 (based on raingauge and radar data) :srcset: /auto_examples/images/sphx_glr_qpe_manila_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyresample/geometry.py:1617: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_params = self.crs.to_proj4() /opt/conda/envs/swirlspy/lib/python3.6/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 36.536 seconds) .. _sphx_glr_download_auto_examples_qpe_manila.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: qpe_manila.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: qpe_manila.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_