.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/ltg_bin_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_ltg_bin_hk.py: Simple Lightning Nowcast (Hong Kong) ======================================================== This example demonstrates how to perform a simple lightning nowcast by the extrapolation of lightning strikes, using data from Hong Kong. .. GENERATED FROM PYTHON SOURCE LINES 10-13 Definitions -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 15-16 Import all required modules and methods: .. GENERATED FROM PYTHON SOURCE LINES 16-66 .. 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 import pyproj 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 # Python com-swirls package to calculate motion field (rover) and semi-lagrangian advection from swirlspy.qpf import rover, sla # 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 lightning function from swirlspy.obs import Lightning from swirlspy.ltg.map import gen_ltg_binary # swirlspy standardize data function from swirlspy.utils import standardize_attr, FrameType # directory constants from swirlspy.tests.samples import DATA_DIR from swirlspy.tests.outputs import OUTPUT_DIR warnings.filterwarnings("ignore") plt.switch_backend('agg') .. GENERATED FROM PYTHON SOURCE LINES 67-72 Initialising ------------------------------------------------------------------------- Define the target grid: .. GENERATED FROM PYTHON SOURCE LINES 72-94 .. code-block:: default # Defining target grid area_id = "hk1980_250km" description = ("A 1km 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 = get_area_def( area_id, description, proj_id, projection, x_size, y_size, area_extent ) .. GENERATED FROM PYTHON SOURCE LINES 95-96 Define basemap: .. GENERATED FROM PYTHON SOURCE LINES 96-130 .. code-block:: default # Load the shape of Hong Kong map_shape_file = os.path.abspath(os.path.join( DATA_DIR, 'shape/hk_hk1980' )) # coastline and province map_with_province = cfeature.ShapelyFeature( list(shpreader.Reader(map_shape_file).geometries()), area_def.to_cartopy_crs() ) # define the plot function def plot_base(ax: plt.Axes, extents: list, crs: ccrs.Projection): # fake the ocean color ax.imshow( np.tile(np.array([[[178, 208, 254]]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=crs, extent=extents, 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 131-133 Defining a basetime, nowcast interval and nowcast duration: .. GENERATED FROM PYTHON SOURCE LINES 133-139 .. code-block:: default basetime = pd.Timestamp('201904201200').floor('min') basetime_str = basetime.strftime('%Y%m%d%H%M') nowcast_interval = pd.Timedelta(30, 'm') nowcast_duration = pd.Timedelta(2, 'h') .. GENERATED FROM PYTHON SOURCE LINES 140-145 Loading Lightning and Radar Data -------------------------------------------------------------------------- Extract Lightning Location Information System (LLIS) files to read. .. GENERATED FROM PYTHON SOURCE LINES 145-168 .. code-block:: default # The interval between LLIS files llis_file_interval = pd.Timedelta(1, 'm') # Finding the string representation # of timestamps marking the relevant LLIS files timestamps = timestamps_ending( basetime, duration=nowcast_interval, interval=llis_file_interval, format='%Y%m%d%H%M' ) # Locating files llis_dir = os.path.join(DATA_DIR, 'llis') located_files = [] for timestamp in timestamps: located_file = locate_file(llis_dir, timestamp) located_files.append(located_file) if located_file is not None: print(f"LLIS file for {timestamp} has been located") .. rst-class:: sphx-glr-script-out .. code-block:: none LLIS file for 201904201130 has been located LLIS file for 201904201131 has been located LLIS file for 201904201132 has been located LLIS file for 201904201133 has been located LLIS file for 201904201134 has been located LLIS file for 201904201135 has been located LLIS file for 201904201136 has been located LLIS file for 201904201137 has been located LLIS file for 201904201138 has been located LLIS file for 201904201139 has been located LLIS file for 201904201140 has been located LLIS file for 201904201141 has been located LLIS file for 201904201142 has been located LLIS file for 201904201143 has been located LLIS file for 201904201144 has been located LLIS file for 201904201145 has been located LLIS file for 201904201146 has been located LLIS file for 201904201147 has been located LLIS file for 201904201148 has been located LLIS file for 201904201149 has been located LLIS file for 201904201150 has been located LLIS file for 201904201151 has been located LLIS file for 201904201152 has been located LLIS file for 201904201153 has been located LLIS file for 201904201154 has been located LLIS file for 201904201155 has been located LLIS file for 201904201156 has been located LLIS file for 201904201157 has been located LLIS file for 201904201158 has been located LLIS file for 201904201159 has been located .. GENERATED FROM PYTHON SOURCE LINES 169-171 Read data from files into Lightning objects. .. GENERATED FROM PYTHON SOURCE LINES 171-181 .. code-block:: default # Initialising Lightning objects # Coordinates are geodetic ltg_object_geodetic = Lightning( located_files, 'geodetic', basetime - nowcast_interval, basetime ) .. GENERATED FROM PYTHON SOURCE LINES 182-185 Locating IRIS REF radar files needed to generate motion field. The IRIS REF files are marked by start time. .. GENERATED FROM PYTHON SOURCE LINES 185-211 .. code-block:: default # The interval between the radar files rad_interval = pd.Timedelta(6, 'm') # Finding the string representation # of the timestamps # marking the required files # In this case, the most recent radar files # separated by `nowcast_interval` is required rad_timestamps = timestamps_ending( basetime, duration=nowcast_interval + rad_interval, interval=rad_interval ) rad_timestamps = [rad_timestamps[0], rad_timestamps[-1]] rad_dir = os.path.join(DATA_DIR, 'iris') located_rad_files = [] for timestamp in rad_timestamps: located_file = locate_file(rad_dir, timestamp) located_rad_files.append(located_file) if located_file is not None: print(f"IRIS REF file for {timestamp} has been located") .. rst-class:: sphx-glr-script-out .. code-block:: none IRIS REF file for 1904201124 has been located IRIS REF file for 1904201154 has been located .. GENERATED FROM PYTHON SOURCE LINES 212-215 Read from iris radar files using read_iris_grid(). The data is in the Azimuthal Equidistant Projection. .. GENERATED FROM PYTHON SOURCE LINES 215-222 .. code-block:: default reflec_unproj_list = [] for file in located_rad_files: reflec = read_iris_grid(file) reflec_unproj_list.append(reflec) .. GENERATED FROM PYTHON SOURCE LINES 223-231 Data Reprojection --------------------------------------------------------------------------- Since the coordinates of the Lightning objects are geodetic, a conversion from geodetic coordinates to that of the user-defined projection system is required. This can be achieved by `swirlspy.obs.Lightning.reproject()` .. GENERATED FROM PYTHON SOURCE LINES 231-240 .. code-block:: default # pyproj representation of the map projection of the input data inProj = pyproj.Proj(init='epsg:4326') # pyproj representation of the intended map projection outProj = pyproj.Proj(area_def.proj_str) # reproject ltg_object = ltg_object_geodetic.reproject(inProj, outProj, 'HK1980') .. GENERATED FROM PYTHON SOURCE LINES 241-246 Radar data also needs to be reprojected from Azimuthal Equidistant to HK1980. This is achieved by the `swirlspy.core.resample` function. Following reprojection, the individual xarray.DataArrays can be concatenated to a single xarray.DataArrays. .. GENERATED FROM PYTHON SOURCE LINES 246-262 .. code-block:: default reflec_list = [] for z in reflec_unproj_list: reflec = grid_resample( z, z.attrs['area_def'], area_def, coord_label=['easting', 'northing'] ) reflec_list.append(reflec) reflectivity = xr.concat(reflec_list, dim='time') standardize_attr(reflectivity, frame_type=FrameType.dBZ) print(reflectivity) .. rst-class:: sphx-glr-script-out .. code-block:: none array([[[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]]]) Coordinates: * time (time) datetime64[ns] 2019-04-20T11:24:00 2019-04-20T11:54:00 * northing (northing) float64 1.068e+06 1.068e+06 ... 5.705e+05 5.695e+05 * easting (easting) float64 5.875e+05 5.885e+05 ... 1.086e+06 1.086e+06 Attributes: long_name: Reflectivity units: dBZ projection: hk1980 site: (114.12333341315389, 22.41166664287448, 580) radar_height: 8 area_def: Area ID: hk1980_250km\nDescription: A 1km resolution recta... proj_site: (830755.3393923986, 830262.1449904075) step_size: 0 days 00:30:00 zero_value: nan frame_type: FrameType.dBZ .. GENERATED FROM PYTHON SOURCE LINES 263-269 Grid lightning observations ----------------------------------------------------------------- Generate a binary grid of lighting observations, where 0 indicates no lightning and 1 indicates that there is lightning. .. GENERATED FROM PYTHON SOURCE LINES 269-287 .. code-block:: default ltg_bin = gen_ltg_binary(ltg_object, area_def, coord_label=['easting', 'northing']) print(ltg_bin) # Duplicating ltg_density along the time dimension # so that it can be advected. ltg_bin1 = ltg_bin.copy() ltg_bin1.coords['time'] = [basetime - nowcast_interval] ltg_bin_concat = xr.concat( [ltg_bin1, ltg_bin], dim='time' ) # Identifying zero value ltg_bin_concat.attrs['zero_value'] = np.nan .. rst-class:: sphx-glr-script-out .. code-block:: none array([[[False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], ..., [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False]]]) Coordinates: * time (time) datetime64[ns] 2019-04-20T12:00:00 * northing (northing) float64 1.068e+06 1.068e+06 ... 5.705e+05 5.695e+05 * easting (easting) float64 5.875e+05 5.885e+05 ... 1.086e+06 1.086e+06 Attributes: long_name: Lightning Occurrence units: 0/1 .. GENERATED FROM PYTHON SOURCE LINES 288-295 Generating motion fields and extrapolating -------------------------------------------------------------- #. Generate motion field using ROVER. #. Extrapolate gridded lightning strikes using the motion field by Semi-Lagrangian Advection. .. GENERATED FROM PYTHON SOURCE LINES 295-304 .. code-block:: default motion = rover(reflectivity) # rover steps = int(nowcast_duration.total_seconds() // nowcast_interval.total_seconds()) ltg_frames = sla(ltg_bin_concat, motion, nowcast_steps=steps-1) # sla print(ltg_frames) .. rst-class:: sphx-glr-script-out .. code-block:: none RUNNING 'rover' FOR EXTRAPOLATION..... array([[[ 0., 0., 0., ..., 0., 0., 0.], [ 0., 0., 0., ..., 0., 0., 0.], [ 0., 0., 0., ..., 0., 0., 0.], ..., [ 0., 0., 0., ..., 0., 0., 0.], [ 0., 0., 0., ..., 0., 0., 0.], [ 0., 0., 0., ..., 0., 0., 0.]], [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]]]) Coordinates: * time (time) datetime64[ns] 2019-04-20T12:00:00 ... 2019-04-20T13:30:00 * northing (northing) float64 1.068e+06 1.068e+06 ... 5.705e+05 5.695e+05 * easting (easting) float64 5.875e+05 5.885e+05 ... 1.086e+06 1.086e+06 Attributes: long_name: Lightning Occurrence units: 0/1 zero_value: nan .. GENERATED FROM PYTHON SOURCE LINES 305-311 Visualization --------------------------------------------------------------- Now, let us visualize the lightning nowcast up to a duration of 2 hours. .. GENERATED FROM PYTHON SOURCE LINES 311-417 .. code-block:: default # Colormap and colorbar parameters n_col = 2 cmap = ListedColormap(['#FFFFFF', '#000000']) levels = [0., 0.5, 1.] ticks = np.array([0, 1]) norm = BoundaryNorm(levels, cmap.N, clip=True) mappable = ScalarMappable(cmap=cmap, norm=norm) mappable.set_array([]) # Defining the crs crs = area_def.to_cartopy_crs() # Defining area x = ltg_frames.coords['easting'].values y = ltg_frames.coords['northing'].values extents = [ np.min(x), np.max(x), np.min(y), np.max(y) ] qx = motion.coords['easting'].values[::20] qy = motion.coords['northing'].values[::20] qu = motion.values[0, ::20, ::20] qv = motion.values[1, ::20, ::20] fig = 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(ltg_frames.time.values): time = pd.Timestamp(t) row = i // 2 col = i % 2 ax = fig.add_subplot(gs[row, col], projection=crs) # plot base map plot_base(ax, extents, crs) # plot lightning frame = ltg_frames.sel(time=time) frame.where(frame > levels[1]).plot( cmap=cmap, norm=norm, add_colorbar=False ) # plot motion vector ax.quiver( qx, qy, qu, qv, pivot='mid' ) ax.set_title('') ax.text( extents[0], extents[3], textwrap.dedent( """ Lightning Density Based @ {baseTime} """ ).format( baseTime=basetime.strftime('%H:%MH') ).strip(), fontsize=10, va='bottom', ha='left', linespacing=1 ) ax.text( extents[1], extents[3], textwrap.dedent( """ {validDate} Valid @ {validTime} """ ).format( validDate=basetime.strftime('%Y-%m-%d'), validTime=time.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) tick_locs = (np.arange(n_col) + 0.5)*(n_col-1)/n_col cbar.set_ticks(tick_locs) cbar.set_ticklabels(np.arange(n_col)) cbar.ax.set_ylabel( f"{ltg_bin.attrs['long_name']}[{ltg_bin.attrs['units']}]" ) fig.savefig( os.path.join(OUTPUT_DIR, "ltg-bin-hk.png"), bbox_inches='tight' ) .. image-sg:: /auto_examples/images/sphx_glr_ltg_bin_hk_001.png :alt: ltg bin hk :srcset: /auto_examples/images/sphx_glr_ltg_bin_hk_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 52.938 seconds) .. _sphx_glr_download_auto_examples_ltg_bin_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: ltg_bin_hk.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ltg_bin_hk.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_