.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/sprog_hk.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_sprog_hk.py: SPROG (Hong Kong) ======================================================== This example demonstrates how to use SPROG to forecast rainfall up to three hours, using rain guage and radar data from Hong Kong. .. GENERATED FROM PYTHON SOURCE LINES 10-13 Setup -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 15-16 Import all required modules and methods: .. GENERATED FROM PYTHON SOURCE LINES 16-40 .. code-block:: default import os import numpy as np import pandas as pd import xarray as xr import textwrap from pyresample import utils import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import cartopy.io.shapereader as shpreader from matplotlib.gridspec import GridSpec from matplotlib.colors import BoundaryNorm, LinearSegmentedColormap from matplotlib.colors import ListedColormap from matplotlib.cm import ScalarMappable from swirlspy.rad.iris import read_iris_grid from swirlspy.qpe.utils import locate_file, timestamps_ending from swirlspy.core.resample import grid_resample from swirlspy.utils import FrameType, standardize_attr, FrameType, conversion from swirlspy.qpf import sprog, dense_lucaskanade .. GENERATED FROM PYTHON SOURCE LINES 41-42 Define working directory and nowcast parameters: .. GENERATED FROM PYTHON SOURCE LINES 42-52 .. code-block:: default # working_dir = os.path.join(os.getcwd(), 'swirlspy/examples') working_dir = os.getcwd() radar_dir = os.path.abspath( os.path.join(working_dir, '../tests/samples/iris/ppi') ) # Set nowcast parameters n_timesteps = int(3 * 60 / 6) # 3 hours, each timestamp is 6 minutes .. GENERATED FROM PYTHON SOURCE LINES 53-54 Define the user grid: .. GENERATED FROM PYTHON SOURCE LINES 54-74 .. code-block:: default 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_tgt = 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/sprog_hk.py:71: 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 75-76 Define the base map: .. GENERATED FROM PYTHON SOURCE LINES 76-110 .. code-block:: default # Load the shape of Hong Kong map_shape_file = os.path.abspath(os.path.join( working_dir, '../tests/samples/shape/hk' )) # coastline and province map_with_province = cfeature.ShapelyFeature( list(shpreader.Reader(map_shape_file).geometries()), ccrs.PlateCarree() ) # define the plot function def plot_base(ax: plt.Axes, extents: list, crs: ccrs.Projection): ax.set_extent(extents, crs=crs) # fake the ocean color ax.imshow( np.tile(np.array([[[178, 208, 254]]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180], zorder=-1 ) # coastline, province and state, color ax.add_feature( map_with_province, facecolor=cfeature.COLORS['land'], edgecolor='none', zorder=0 ) # overlay coastline, province and state without color ax.add_feature( map_with_province, facecolor='none', edgecolor='gray', linewidth=0.5 ) ax.set_title('') .. GENERATED FROM PYTHON SOURCE LINES 111-112 Log the start time for reference: .. GENERATED FROM PYTHON SOURCE LINES 112-117 .. code-block:: default start_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 118-121 Loading Radar Data --------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 121-173 .. code-block:: default # Specify the basetime basetime = pd.Timestamp('201902190800') # Generate a list of timestamps of the radar data files located_files = [] radar_ts = timestamps_ending( basetime, duration=pd.Timedelta(minutes=60), exclude_end=False ) for timestamp in radar_ts: located_files.append(locate_file(radar_dir, timestamp)) # Read in the radar data reflectivity_list = [] # stores reflec from read_iris_grid() for filename in located_files: reflec = read_iris_grid(filename) reflectivity_list.append(reflec) # Reproject the radar data to the user-defined grid area_def_src = reflectivity_list[0].attrs['area_def'] reproj_reflectivity_list = [] for reflec in reflectivity_list: reproj_reflec = grid_resample( reflec, area_def_src, area_def_tgt, coord_label=['x', 'y'] ) reproj_reflectivity_list.append(reproj_reflec) # Standardize reflectivity xarrays raw_frames = xr.concat(reproj_reflectivity_list, dim='time').sortby(['y'], ascending=False) standardize_attr(raw_frames, frame_type=FrameType.dBZ) # Transform from reflecitiy to rainfall rate frames = conversion.to_rainfall_rate(raw_frames, True, a=58.53, b=1.56) # Set the fill value frames.attrs['zero_value'] = -15.0 # Apply threshold to -10dBR i.e. 0.1mm/h threshold = -10.0 frames.values[frames.values < threshold] = frames.attrs['zero_value'] # Set missing values with the fill value frames.values[~np.isfinite(frames.values)] = frames.attrs['zero_value'] # Log the time for record initialising_time = pd.Timestamp.now() .. rst-class:: sphx-glr-script-out .. code-block:: pytb Traceback (most recent call last): File "/tmp/build/docs/swirlspy/swirlspy/examples/sprog_hk.py", line 138, in reflec = read_iris_grid(filename) File "/tmp/build/docs/swirlspy/swirlspy/rad/_iris.py", line 407, in read_iris_grid raise ValueError("Invalid file") from e ValueError: Invalid file .. GENERATED FROM PYTHON SOURCE LINES 174-177 Running Lucas Kanade Optical flow and S-PROG -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 177-197 .. code-block:: default # Estimate the motion field motion = dense_lucaskanade(frames) motion_time = pd.Timestamp.now() # Generate forecast rainrate field forcast_frames = sprog( frames, motion, n_timesteps, n_cascade_levels=8, R_thr=threshold, decomp_method="fft", bandpass_filter_method="gaussian", probmatching_method="mean", ) sprog_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 198-203 Generating radar reflectivity maps ---------------------------------- Define the color scale and format of the plot. .. GENERATED FROM PYTHON SOURCE LINES 203-317 .. code-block:: default # Defining colour scale and format. levels = [ -32768, 10, 15, 20, 24, 28, 32, 34, 38, 41, 44, 47, 50, 53, 56, 58, 60, 62 ] cmap = ListedColormap([ '#FFFFFF00', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433', '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200', '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0' ]) norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) mappable = ScalarMappable(cmap=cmap, norm=norm) mappable.set_array([]) # Defining the crs crs = area_def_tgt.to_cartopy_crs() # Defining area x = frames.coords['x'].values y = frames.coords['y'].values x_d = x[1] - x[0] y_d = y[1] - y[0] extents = [x[0], y[0], x[-1], y[-1]] # Generating a time steps for every hour time_steps = [ basetime + pd.Timedelta(minutes=6*i) for i in range(n_timesteps + 1) if i % 10 == 0 ] ref_frames = conversion.to_reflectivity(forcast_frames, True) ref_frames.data[ref_frames.data < 0.1] = np.nan ref_frames = xr.concat([raw_frames[:-1, ...], ref_frames], dim='time') ref_frames.attrs['values_name'] = 'Reflectivity 2km CAPPI' standardize_attr(ref_frames) qx = motion.coords['x'].values[::5] qy = motion.coords['y'].values[::5] qu = motion.values[0, ::5, ::5] qv = motion.values[1, ::5, ::5] fig: plt.Figure = plt.figure(figsize=(8, 8), frameon=False) gs = GridSpec( 2, 2, figure=fig, wspace=0.03, hspace=-0.25, top=0.95, bottom=0.05, left=0.17, right=0.845 ) for i, t in enumerate(time_steps): row = i // 2 col = i % 2 ax = fig.add_subplot(gs[row, col], projection=crs) # plot base map plot_base(ax, extents, crs) # plot reflectivity frame = ref_frames.sel(time=t) im = ax.imshow(frame.values, cmap=cmap, norm=norm, interpolation='nearest', extent=extents) # plot motion vector ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20) ax.text( extents[0], extents[1], textwrap.dedent( """ Reflectivity Based @ {baseTime} """ ).format( baseTime=basetime.strftime('%H:%MH') ).strip(), fontsize=10, va='bottom', ha='left', linespacing=1 ) ax.text( extents[2] - (extents[2] - extents[0]) * 0.03, extents[1], textwrap.dedent( """ {validDate} Valid @ {validTime} """ ).format( validDate=basetime.strftime('%Y-%m-%d'), validTime=t.strftime('%H:%MH') ).strip(), fontsize=10, va='bottom', ha='right', linespacing=1 ) cbar_ax = fig.add_axes([0.875, 0.125, 0.03, 0.75]) cbar = fig.colorbar( mappable, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g') cbar.ax.set_ylabel(ref_frames.attrs['values_name'], rotation=90) fig.savefig( os.path.join( working_dir, "../tests/outputs/sprog-reflectivity.png" ), bbox_inches='tight' ) radar_image_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 318-324 Accumulating hourly rainfall for 3-hour forecast ------------------------------------------------ Hourly accumulated rainfall is calculated every 30 minutes, the first endtime is the basetime i.e. T+30min. .. GENERATED FROM PYTHON SOURCE LINES 324-343 .. code-block:: default # Optional, convert to rainfall depth rf_frames = conversion.to_rainfall_depth(ref_frames, a=58.53, b=1.56) # Compute hourly accumulated rainfall every 60 minutes. acc_rf_frames = conversion.acc_rainfall_depth( rf_frames, basetime, basetime + pd.Timedelta(hours=3), pd.Timedelta(minutes=60) ) # Replace zero value with NaN acc_rf_frames.data[acc_rf_frames.data <= acc_rf_frames.attrs['zero_value']] = np.nan acc_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 344-347 Generating radar reflectivity maps ----------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 347-434 .. code-block:: default # Defining colour scale and format. levels = [ 0, 0.5, 2, 5, 10, 20, 30, 40, 50, 70, 100, 150, 200, 300, 400, 500, 600, 700 ] cmap = ListedColormap([ '#ffffff00', '#9bf7f7', '#00ffff', '#00d5cc', '#00bd3d', '#2fd646', '#9de843', '#ffdd41', '#ffac33', '#ff621e', '#d23211', '#9d0063', '#e300ae', '#ff00ce', '#ff57da', '#ff8de6', '#ffe4fd' ]) norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) mappable = ScalarMappable(cmap=cmap, norm=norm) mappable.set_array([]) fig: plt.Figure = plt.figure(figsize=(8, 8), frameon=False) gs = GridSpec( 2, 2, figure=fig, wspace=0.03, hspace=-0.25, top=0.95, bottom=0.05, left=0.17, right=0.845 ) for i, t in enumerate(acc_rf_frames.coords['time'].values): row = i // 2 col = i % 2 ax = fig.add_subplot(gs[row, col], projection=crs) # plot base map plot_base(ax, extents, crs) # plot accumulated rainfall depth t = pd.Timestamp(t) frame = acc_rf_frames.sel(time=t) im = ax.imshow(frame.values, cmap=cmap, norm=norm, interpolation='nearest', extent=extents) ax.text( extents[0], extents[1], textwrap.dedent( """ Hourly Rainfall Based @ {baseTime} """ ).format( baseTime=basetime.strftime('%H:%MH') ).strip(), fontsize=10, va='bottom', ha='left', linespacing=1 ) ax.text( extents[2] - (extents[2] - extents[0]) * 0.03, extents[1], textwrap.dedent( """ {validDate} Valid @ {validTime} """ ).format( validDate=basetime.strftime('%Y-%m-%d'), validTime=t.strftime('%H:%MH') ).strip(), fontsize=10, va='bottom', ha='right', linespacing=1 ) cbar_ax = fig.add_axes([0.875, 0.125, 0.03, 0.75]) cbar = fig.colorbar( mappable, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g') cbar.ax.set_ylabel(acc_rf_frames.attrs['values_name'], rotation=90) fig.savefig( os.path.join( working_dir, "../tests/outputs/sprog-rainfall.png" ), bbox_inches='tight' ) ptime = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 435-438 Checking run time of each component -------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 438-455 .. code-block:: default print(f"Start time: {start_time}") print(f"Initialising time: {initialising_time}") print(f"Motion field time: {motion_time}") print(f"S-PROG time: {sprog_time}") print(f"Plotting radar image time: {radar_image_time}") print(f"Accumulating rainfall time: {acc_time}") print(f"Plotting rainfall maps: {ptime}") print(f"Time to initialise: {initialising_time - start_time}") print(f"Time to run motion field: {motion_time - initialising_time}") print(f"Time to perform S-PROG: {sprog_time - motion_time}") print(f"Time to plot radar image: {radar_image_time - sprog_time}") print(f"Time to accumulate rainfall: {acc_time - radar_image_time}") print(f"Time to plot rainfall maps: {ptime - acc_time}") print(f"Total: {ptime - start_time}") .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.061 seconds) .. _sphx_glr_download_auto_examples_sprog_hk.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: sprog_hk.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: sprog_hk.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_