.. 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-12 Definitions ----------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 12-42 .. code-block:: default import os from copy import copy import textwrap import pyproj import numpy as np import pandas as pd import xarray as xr import matplotlib.pyplot as plt from matplotlib.colors import BoundaryNorm, ListedColormap from matplotlib.gridspec import GridSpec import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.io import shapereader from matplotlib.cm import ScalarMappable from pyresample import utils from swirlspy.qpf import rover from swirlspy.qpf import sla from swirlspy.qpe.utils import timestamps_ending, locate_file from swirlspy.rad.iris import read_iris_grid from swirlspy.core.resample import grid_resample from swirlspy.obs import Lightning from swirlspy.ltg.map import gen_ltg_binary from swirlspy.utils import standardize_attr, FrameType plt.switch_backend('agg') THIS_DIR = os.getcwd() .. GENERATED FROM PYTHON SOURCE LINES 43-48 Initialising ------------------------------------------------------------------------- Define the target grid: .. GENERATED FROM PYTHON SOURCE LINES 48-70 .. 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 = 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/ltg_bin_hk.py:67: 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 71-72 Define basemap: .. GENERATED FROM PYTHON SOURCE LINES 72-106 .. code-block:: default # Load the shape of Hong Kong map_shape_file = os.path.abspath(os.path.join( THIS_DIR, '../tests/samples/shape/hk_hk1980' )) # coastline and province map_with_province = cfeature.ShapelyFeature( list(shapereader.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('') .. 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() .. GENERATED FROM PYTHON SOURCE LINES 107-109 Defining a basetime, nowcast interval and nowcast duration: .. GENERATED FROM PYTHON SOURCE LINES 109-116 .. 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 117-122 Loading Lightning and Radar Data -------------------------------------------------------------------------- Extract Lightning Location Information System (LLIS) files to read. .. GENERATED FROM PYTHON SOURCE LINES 122-145 .. 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(THIS_DIR, '../tests/samples/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 146-148 Read data from files into Lightning objects. .. GENERATED FROM PYTHON SOURCE LINES 148-158 .. 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 159-162 Locating IRIS REF radar files needed to generate motion field. The IRIS REF files are marked by start time. .. GENERATED FROM PYTHON SOURCE LINES 162-188 .. 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(THIS_DIR, '../tests/samples/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 189-192 Read from iris radar files using read_iris_grid(). The data is in the Azimuthal Equidistant Projection. .. GENERATED FROM PYTHON SOURCE LINES 192-199 .. code-block:: default reflec_unproj_list = [] for file in located_rad_files: reflec = read_iris_grid(file) reflec_unproj_list.append(reflec) .. rst-class:: sphx-glr-script-out .. code-block:: pytb Traceback (most recent call last): File "/tmp/build/docs/swirlspy/swirlspy/examples/ltg_bin_hk.py", line 195, in reflec = read_iris_grid(file) 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 200-208 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 208-217 .. 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 218-223 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 223-239 .. 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) .. GENERATED FROM PYTHON SOURCE LINES 240-246 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 246-264 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 265-272 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 272-281 .. 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) .. GENERATED FROM PYTHON SOURCE LINES 282-288 Visualization --------------------------------------------------------------- Now, let us visualize the lightning nowcast up to a duration of 2 hours. .. GENERATED FROM PYTHON SOURCE LINES 288-396 .. 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( THIS_DIR, "../tests/outputs/ltg-bin-hk.png" ), bbox_inches='tight' ) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.802 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 `_