.. 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-12 Definitions -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 14-15 Import all required modules and methods: .. GENERATED FROM PYTHON SOURCE LINES 15-60 .. code-block:: default # Python package to allow system command line functions import os # Python package to manage warning message import warnings # Python package for time calculations import pandas as pd # Python package for numerical calculations import numpy as np # Python package for xarrays to read and handle netcdf data import xarray as xr # Python package for text formatting import textwrap # Python package for projection description from pyresample import get_area_def # Python package for projection import cartopy.crs as ccrs # Python package for land/sea features import cartopy.feature as cfeature # Python package for reading map shape file import cartopy.io.shapereader as shpreader # Python package for creating plots from matplotlib import pyplot as plt # Python package for output grid format from matplotlib.gridspec import GridSpec # Python package for colorbars from matplotlib.colors import BoundaryNorm, ListedColormap from matplotlib.cm import ScalarMappable # swirlspy data parser function from swirlspy.rad.iris import read_iris_grid # swirlspy test data source locat utils function from swirlspy.qpe.utils import timestamps_ending, locate_file # swirlspy regrid function from swirlspy.core.resample import grid_resample # swirlspy standardize data function from swirlspy.utils import FrameType, standardize_attr, conversion # swirlspy pysteps integrated package from swirlspy.qpf import steps, dense_lucaskanade # directory constants from swirlspy.tests.samples import DATA_DIR from swirlspy.tests.outputs import OUTPUT_DIR warnings.filterwarnings("ignore") .. GENERATED FROM PYTHON SOURCE LINES 61-62 Define the working directory and nowcast parameters .. GENERATED FROM PYTHON SOURCE LINES 62-71 .. code-block:: default radar_dir = os.path.abspath( os.path.join(DATA_DIR, '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 72-73 Define the User Grid .. GENERATED FROM PYTHON SOURCE LINES 73-91 .. 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 = get_area_def( area_id, description, proj_id, projection, x_size, y_size, area_extent ) .. GENERATED FROM PYTHON SOURCE LINES 92-93 Define the plotting function: .. GENERATED FROM PYTHON SOURCE LINES 93-126 .. code-block:: default # Defining plot parameters map_shape_file = os.path.abspath(os.path.join( DATA_DIR, '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 127-129 Loading radar data --------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 129-185 .. 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() .. GENERATED FROM PYTHON SOURCE LINES 186-189 Running Lucas Kanade Optical flow and S-PROG ------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 189-216 .. 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() .. rst-class:: sphx-glr-script-out .. code-block:: none Computing STEPS nowcast: ------------------------ Inputs: ------- input dimensions: 500x500 km/pixel: 2 time step: 6 minutes Methods: -------- extrapolation: semilagrangian bandpass filter: gaussian decomposition: fft noise generator: nonparametric noise adjustment: no velocity perturbator: bps conditional statistics: no precip. mask method: incremental probability matching: mean FFT method: numpy domain: spatial Parameters: ----------- number of time steps: 30 ensemble size: 3 parallel threads: 1 number of cascade levels: 8 order of the AR(p) model: 2 velocity perturbations, parallel: 10.88,0.23,-7.68 velocity perturbations, perpendicular: 5.76,0.31,-2.72 precip. intensity threshold: -10 ************************************************ * Correlation coefficients for cascade levels: * ************************************************ ----------------------------------------- | Level | Lag-1 | Lag-2 | ----------------------------------------- | 1 | 0.998995 | 0.996974 | ----------------------------------------- | 2 | 0.998526 | 0.996270 | ----------------------------------------- | 3 | 0.993048 | 0.983181 | ----------------------------------------- | 4 | 0.972410 | 0.927100 | ----------------------------------------- | 5 | 0.868747 | 0.695263 | ----------------------------------------- | 6 | 0.515560 | 0.239395 | ----------------------------------------- | 7 | 0.103358 | 0.020485 | ----------------------------------------- | 8 | 0.003543 | 0.011290 | ----------------------------------------- **************************************** * AR(p) parameters for cascade levels: * **************************************** ------------------------------------------------------ | Level | Phi-1 | Phi-2 | Phi-0 | ------------------------------------------------------ | 1 | 1.505001 | -0.506515 | 0.038644 | ------------------------------------------------------ | 2 | 1.264189 | -0.266055 | 0.052323 | ------------------------------------------------------ | 3 | 1.205369 | -0.213808 | 0.114992 | ------------------------------------------------------ | 4 | 1.302659 | -0.339619 | 0.219412 | ------------------------------------------------------ | 5 | 1.079346 | -0.242416 | 0.480483 | ------------------------------------------------------ | 6 | 0.534103 | -0.035967 | 0.856299 | ------------------------------------------------------ | 7 | 0.102334 | 0.009908 | 0.994595 | ------------------------------------------------------ | 8 | 0.003503 | 0.011277 | 0.999930 | ------------------------------------------------------ Starting nowcast computation. Computing nowcast for time step 1... done. Computing nowcast for time step 2... done. Computing nowcast for time step 3... done. Computing nowcast for time step 4... done. Computing nowcast for time step 5... done. Computing nowcast for time step 6... done. Computing nowcast for time step 7... done. Computing nowcast for time step 8... done. Computing nowcast for time step 9... done. Computing nowcast for time step 10... done. Computing nowcast for time step 11... done. Computing nowcast for time step 12... done. Computing nowcast for time step 13... done. Computing nowcast for time step 14... done. Computing nowcast for time step 15... done. Computing nowcast for time step 16... done. Computing nowcast for time step 17... done. Computing nowcast for time step 18... done. Computing nowcast for time step 19... done. Computing nowcast for time step 20... done. Computing nowcast for time step 21... done. Computing nowcast for time step 22... done. Computing nowcast for time step 23... done. Computing nowcast for time step 24... done. Computing nowcast for time step 25... done. Computing nowcast for time step 26... done. Computing nowcast for time step 27... done. Computing nowcast for time step 28... done. Computing nowcast for time step 29... done. Computing nowcast for time step 30... done. .. GENERATED FROM PYTHON SOURCE LINES 217-220 Generating radar reflectivity maps ----------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 220-341 .. 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( OUTPUT_DIR, "steps-reflectivity.png" ), bbox_inches='tight' ) radar_image_time = pd.Timestamp.now() .. image-sg:: /auto_examples/images/sphx_glr_steps_hk_001.png :alt: steps hk :srcset: /auto_examples/images/sphx_glr_steps_hk_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 342-348 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 348-369 .. 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 370-373 Generating radar reflectivity maps ----------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 373-465 .. 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( OUTPUT_DIR, "steps-rainfall.png" ), bbox_inches='tight' ) ptime = pd.Timestamp.now() .. image-sg:: /auto_examples/images/sphx_glr_steps_hk_002.png :alt: steps hk :srcset: /auto_examples/images/sphx_glr_steps_hk_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 466-469 Checking run time of each component -------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 469-486 .. 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-script-out .. code-block:: none Start time: 2024-04-18 03:05:50.706837 Initialising time: 2024-04-18 03:05:52.834310 Motion field time: 2024-04-18 03:05:54.147495 STEPS time: 2024-04-18 03:06:22.945160 Plotting radar image time: 2024-04-18 03:06:58.087656 Accumulating rainfall time: 2024-04-18 03:07:00.629242 Plotting rainfall maps: 2024-04-18 03:07:43.447035 Time to initialise: 0 days 00:00:02.127473 Time to run motion field: 0 days 00:00:01.313185 Time to perform STEPS: 0 days 00:00:28.797665 Time to plot radar image: 0 days 00:00:35.142496 Time to accumulate rainfall: 0 days 00:00:02.541586 Time to plot rainfall maps: 0 days 00:00:42.817793 Total: 0 days 00:01:52.740198 .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 1 minutes 54.332 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 `_