.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/steps_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_steps_hk.py: STEPS (Hong Kong) ======================================================== This example demonstrates how to make use of STEPS to perform 3-hour rainfall nowcasting with radar data in Hong Kong .. GENERATED FROM PYTHON SOURCE LINES 9-11 Setup -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 13-14 Import all required modules and methods: .. GENERATED FROM PYTHON SOURCE LINES 14-41 .. code-block:: default import os import textwrap import numpy as np import pandas as pd import xarray as xr 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, 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 from swirlspy.utils import standardize_attr, FrameType from swirlspy.utils import conversion from swirlspy.qpf import steps from swirlspy.qpf import dense_lucaskanade .. GENERATED FROM PYTHON SOURCE LINES 42-43 Define the working directory and nowcast parameters .. GENERATED FROM PYTHON SOURCE LINES 43-54 .. 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 n_ens_members = 3 .. GENERATED FROM PYTHON SOURCE LINES 55-56 Define the User Grid .. GENERATED FROM PYTHON SOURCE LINES 56-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/steps_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 plotting function: .. GENERATED FROM PYTHON SOURCE LINES 76-109 .. code-block:: default # Defining plot parameters 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() ) 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 110-112 Loading radar data --------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 112-168 .. code-block:: default # Log the start time start_time = pd.Timestamp.now() # Define the basetime basetime = pd.Timestamp('201902190800') # Generate the timestamps and locate the files located_files = [] radar_ts = timestamps_ending( basetime, duration=pd.Timedelta(minutes=12), exclude_end=False ) for timestamp in radar_ts: located_files.append(locate_file(radar_dir, timestamp)) # Read in the 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) # Fill in all fields of the xarray of reflectivity data frames = xr.concat(reproj_reflectivity_list, dim='time').sortby(['y'], ascending=False) standardize_attr(frames, frame_type=FrameType.dBZ) # Convert from reflectivity to rainfall rain frames = conversion.to_rainfall_rate(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/steps_hk.py", line 133, 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 169-172 Running Lucas Kanade Optical flow and S-PROG ------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 172-199 .. code-block:: default # Estimate the motion field with Lucas Kanade motion = dense_lucaskanade(frames) # Log the time for record motion_time = pd.Timestamp.now() # Nowcast using STEP forcast_frames = steps( frames, motion, n_timesteps, n_ens_members=n_ens_members, n_cascade_levels=8, R_thr=threshold, kmperpixel=2.0, decomp_method="fft", bandpass_filter_method="gaussian", noise_method="nonparametric", probmatching_method="mean", vel_pert_method="bps", mask_method="incremental", seed=24 ) steps_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 200-203 Generating radar reflectivity maps ----------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 203-324 .. code-block:: default # Defining the colour scale 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 extents = [ x[0], y[0], x[-1], y[-1] ] # Generate 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 qx = motion.coords['x'].values[::25] qy = motion.coords['y'].values[::25] qu = motion.values[0, ::25, ::25] qv = motion.values[1, ::25, ::25] n_rows = len(time_steps) fig: plt.Figure = plt.figure( figsize=(n_ens_members * 3.5 + 1, n_rows * 3.5), frameon=False) gs = GridSpec(n_rows, n_ens_members, figure=fig, wspace=0.03, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845) for row in range(n_rows): for col in range(n_ens_members): ax: plt.Axes = fig.add_subplot(gs[row, col], projection=crs) ensemble = ref_frames.coords['ensembles'].values[col] t = time_steps[row] # plot base map plot_base(ax, extents, crs) # plot reflectivity member = ref_frames.sel(ensembles=ensemble) frame = member.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' ) 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.075, 0.03, 0.85]) 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/steps-reflectivity.png" ), bbox_inches='tight' ) radar_image_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 325-331 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 331-352 .. code-block:: default # Convert from rainfall rate to rainfall rf_frames = conversion.to_rainfall_depth(ref_frames, a=58.53, b=1.56) # Compute hourly accumulated rainfall every 30 minutes. acc_rf_frames = [] for ens in rf_frames.coords['ensembles']: af = conversion.acc_rainfall_depth( rf_frames.sel(ensembles=ens).drop('ensembles'), basetime + pd.Timedelta(minutes=60), basetime + pd.Timedelta(hours=3) ) af_ensembles = af.assign_coords(ensembles=ens) acc_rf_frames.append(af_ensembles.expand_dims('ensembles')) acc_rf_frames = xr.concat(acc_rf_frames, dim='ensembles') # 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 353-356 Generating radar reflectivity maps ----------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 356-448 .. code-block:: default # Defining colour scale 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([]) time_steps = acc_rf_frames.coords['time'].values n_rows = len(time_steps) fig: plt.Figure = plt.figure( figsize=(n_ens_members * 4 + 1, n_rows * 4), frameon=False) gs = GridSpec(n_rows, n_ens_members, figure=fig, wspace=0.03, hspace=-0.25, top=0.95, bottom=0.05, left=0.17, right=0.845) for row in range(n_rows): for col in range(n_ens_members): ax = fig.add_subplot(gs[row, col], projection=crs) ensemble = acc_rf_frames.coords['ensembles'].values[col] t = time_steps[row] # plot base map plot_base(ax, extents, crs) # plot accumulated rainfall depth t = pd.Timestamp(t) member = acc_rf_frames.sel(ensembles=ensemble) frame = member.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.095, 0.03, 0.8]) 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/steps-rainfall.png" ), bbox_inches='tight' ) ptime = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 449-452 Checking run time of each component -------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 452-469 .. code-block:: default print(f"Start time: {start_time}") print(f"Initialising time: {initialising_time}") print(f"Motion field time: {motion_time}") print(f"STEPS time: {steps_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 STEPS: {steps_time - motion_time}") print(f"Time to plot radar image: {radar_image_time - steps_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.059 seconds) .. _sphx_glr_download_auto_examples_steps_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: steps_hk.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: steps_hk.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_