.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/ltg_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_hk.py: Lightning Nowcast (Hong Kong) ======================================================== This example demonstrates how to perform lightning nowcasts of up to three hours, using data from Hong Kong. .. GENERATED FROM PYTHON SOURCE LINES 11-13 Definitions ----------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 13-39 .. code-block:: default import os import xarray import time import pandas as pd import matplotlib import imageio import pyproj import matplotlib.pyplot as plt import cartopy.feature as cfeature import cartopy.crs as ccrs from pyresample import utils from copy import copy from matplotlib.colors import LinearSegmentedColormap, BoundaryNorm from swirlspy.qpf.rover import rover from swirlspy.qpf.sla import sla from swirlspy.qpe.utils import timestamps_ending, locate_file from swirlspy.rad.iris import read_iris from swirlspy.obs.lightning import Lightning from swirlspy.ltg.map import gen_ltgv from swirlspy.core.resample import grid_resample THIS_DIR = os.getcwd() .. rst-class:: sphx-glr-script-out .. code-block:: pytb Traceback (most recent call last): File "/tmp/build/docs/swirlspy/swirlspy/examples/ltg_hk.py", line 28, in from swirlspy.qpf.rover import rover File "/tmp/build/docs/swirlspy/swirlspy/qpf/rover.py", line 6, in from rover.rover import rover as rover_api ImportError: libopencv_core.so.3.4: cannot open shared object file: No such file or directory .. GENERATED FROM PYTHON SOURCE LINES 40-45 Initialising ------------------------------------------------------------------------- Step 1: Define the target grid as a pyresample AreaDefinition. .. GENERATED FROM PYTHON SOURCE LINES 45-67 .. code-block:: default # Defining target grid 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 = utils.get_area_def( area_id, description, proj_id, projection, x_size, y_size, area_extent ) .. GENERATED FROM PYTHON SOURCE LINES 68-70 Step 2: Defining a basetime, nowcast interval and nowcast duration. .. GENERATED FROM PYTHON SOURCE LINES 70-77 .. code-block:: default basetime = pd.Timestamp('20180811073000').floor('min') basetime_str = basetime.strftime('%Y%m%d%H%M') interval = pd.Timedelta(minutes=15) nowcast_duration = pd.Timedelta(hours=3) .. GENERATED FROM PYTHON SOURCE LINES 78-81 Step 3: Extracting Lightning Location Information System (LLIS) lightning files to read. .. GENERATED FROM PYTHON SOURCE LINES 81-104 .. code-block:: default # Defining the times of files and the accumulation endtimes of files endtimes = [basetime - interval, basetime] timestamps = [] for time in endtimes: ts = timestamps_ending( time, duration=interval, interval=pd.Timedelta(minutes=1) ) timestamps.append(ts) # Locating files llis_dir = THIS_DIR + '/../tests/samples/llis' located_files = [] for tss in timestamps: lf = [] for timestamp in tss: lf.append(locate_file(llis_dir, timestamp)) located_files.append(lf) .. GENERATED FROM PYTHON SOURCE LINES 105-107 Step 4: Read data from files into Lightning objects. .. GENERATED FROM PYTHON SOURCE LINES 107-117 .. code-block:: default # Initialising Lightning objects # Coordinates are geodetic ltg_object_list_wgs84 = [] for idx, lst in enumerate(located_files): ltg_object_wgs84 = Lightning( lst, 'WGS84', endtimes[idx]-interval, endtimes[idx] ) ltg_object_list_wgs84.append(ltg_object_wgs84) .. GENERATED FROM PYTHON SOURCE LINES 118-119 Step 5: Locating IRIS RAW radar files needed to generate motion field. .. GENERATED FROM PYTHON SOURCE LINES 119-130 .. code-block:: default # Generating timestrings for required times endtimes_str = [] for time in endtimes: endtimes_str.append(time.strftime('%y%m%d%H%M')) rad_dir = THIS_DIR + '/../tests/samples/iris/' located_rad_files = [] for time_str in endtimes_str: located_rad_files.append(locate_file(rad_dir, time_str)) .. GENERATED FROM PYTHON SOURCE LINES 131-135 Step 6: Read from iris radar files using read_iris(). Projection system is specified in the read_iris() function using the AreaDefinition defined above. The resulting xarray.DataArrays are then concatenated along the time axis. .. GENERATED FROM PYTHON SOURCE LINES 135-146 .. code-block:: default reflec_list = [] for file in located_rad_files: reflec = read_iris( file, area_def=area_def, coord_label=['easting', 'northing'] ) reflec_list.append(reflec) reflectivity = xarray.concat(reflec_list, dim='time') .. GENERATED FROM PYTHON SOURCE LINES 147-153 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. .. GENERATED FROM PYTHON SOURCE LINES 153-205 .. code-block:: default # Define a function for neater conversion def convert_ltg( ltg_object_wgs84, area_def ): """ Convert coordinates in object from WGS84 to HK1980. Parameters ----------- ltg_object_wgs84: swirlspy.obs.lightning.Lightning Lightning object with WGS84 coordinates. area_def: pyresample.geometry.AreaDefinition AreaDefinition of the target grid. Returns ------- ltg_object: swirlspy.obs.lightning.Lightning Rain object in HK1980. """ # Defining Proj objects geod = pyproj.Proj(init="epsg:4326") outproj = pyproj.Proj(area_def.proj_str) lx = [] ly = [] data = [] for tup in ltg_object_wgs84.data: lat = tup[1] lon = tup[2] easting, northing = pyproj.transform( geod, outproj, lon, lat ) lx.append(easting) ly.append(northing) data.append((tup[0], northing, easting)) ltg_object = copy(ltg_object_wgs84) ltg_object.lx = lx ltg_object.ly = ly ltg_object.data = data ltg_object.proj = 'HK80' return ltg_object .. GENERATED FROM PYTHON SOURCE LINES 206-215 .. code-block:: default # Calling the function to convert projection system of coordinates ltg_object_list = [] for ltg_object_wgs84 in ltg_object_list_wgs84: ltg_object = convert_ltg(ltg_object_wgs84, area_def) ltg_object_list.append(ltg_object) .. GENERATED FROM PYTHON SOURCE LINES 216-224 Generating lightning potential field ------------------------------------------------------------------ Generate lightning potential fields using gen_ltgv(). Lightning potential decrease normally according to the provided sigma from the site. Lightning potentials are then concatenated along the time axis. .. GENERATED FROM PYTHON SOURCE LINES 224-242 .. code-block:: default # Obtaining lightning potential field from lightning objects ltgv_list = [] for ltg_object in ltg_object_list: ltgvi = gen_ltgv( ltg_object, area_def, ltg_object.start_time, ltg_object.end_time, sigma=20000, coord_label=['easting', 'northing'] ) ltgv_list.append(ltgvi) # Concatenating lightning potential along time dimension ltgv = xarray.concat(ltgv_list, dim='time') .. GENERATED FROM PYTHON SOURCE LINES 243-250 Generating motion fields and extrapolating -------------------------------------------------------------- 1. Generate motion field using ROVER. 2. Extrapolate lightning potential field using the motion field by Semi-Lagrangian Advection .. GENERATED FROM PYTHON SOURCE LINES 250-258 .. code-block:: default # ROVER u, v = rover(reflectivity) # SLA steps = int(nowcast_duration.total_seconds()) // int(interval.total_seconds()) ltgvf = sla(ltgv, u, v, steps=steps) .. GENERATED FROM PYTHON SOURCE LINES 259-268 Plotting results ---------------------------------------------------------------- Plot motion fields with reflectivities and lightning potential nowcasts. Functions can be written to make plotting neater. Step 1: Write function to plot reflectivity with motion fields. .. GENERATED FROM PYTHON SOURCE LINES 268-358 .. code-block:: default def plot_reflectivity( reflectivity, crs, qx, qy, qu, qv, coastline ): """ Custom function for plotting reflectivities with motion fields. Parameters ------------------------------------------ reflectivity: xarray.DataArray Contains data to be plotted crs: cartopy.crs.CRS Defines the projection system. qx: numpy.2darray Contains the x-coordinate values to be plotted. qy: numpy.2darray Contains the y-coordinate values to be plotted qu: numpy.2darray Contains horizontal motion field values. qv: numpy.2darray Contains vertical motion field values. coastline: cartopy.feature Feature (coastline) to be added. """ # Defining colour levels and format levels = [ -32768, 10, 15, 20, 24, 28, 32, 34, 38, 41, 44, 47, 50, 53, 56, 58, 60, 62 ] cmap = matplotlib.colors.ListedColormap([ '#FFFFFF', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433', '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200', '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0' ]) norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) images = [] for time in reflectivity.coords['time'].values: # Plotting plt.figure(figsize=(28, 21)) ax = plt.axes(projection=crs) ax.add_feature(hires) quadmesh = reflectivity.sel(time=time).plot( cmap=cmap, norm=norm, extend='neither' ) cbar = quadmesh.colorbar cbar.ax.set_ylabel( reflectivity.attrs['long_name']+'[' + reflectivity.attrs['units']+']', fontsize=28 ) cbar.ax.tick_params(labelsize=24) ax.quiver( qx, qy, qu, qv, pivot='mid') ax.xaxis.set_visible(True) ax.yaxis.set_visible(True) ax.xaxis.set_tick_params(labelsize=24) ax.yaxis.set_tick_params(labelsize=24) ax.xaxis.label.set_size(28) ax.yaxis.label.set_size(28) pdtime = pd.Timestamp(time) plt.title( "Reflectivity Field with Motion Vectors\n" f"Time: {pdtime}", fontsize=25 ) plt.savefig( THIS_DIR + f"/../tests/outputs/z_{pdtime.strftime('%Y%m%d%H%M')}.png" ) images.append( imageio.imread( THIS_DIR + f"/../tests/outputs/z_{pdtime.strftime('%Y%m%d%H%M')}.png" ) ) imageio.mimsave( THIS_DIR + "/../tests/outputs/reflectivity.gif", images, duration=1 ) .. GENERATED FROM PYTHON SOURCE LINES 359-360 Step 2: Write a function to plot nowcasted lightning potential. .. GENERATED FROM PYTHON SOURCE LINES 360-459 .. code-block:: default # Plotting potential def plot_nowcast( ltgvf, crs, qx, qy, qu, qv, basetime, interval, coastline ): """ Custom function for plotting lightning potential nowcasts with motion fields. Parameters ------------------------------------------ ltgvf: xarray.DataArray Contains data to be plotted crs: cartopy.crs.CRS Defines the projection system. qx: numpy.2darray Contains the x-coordinate values to be plotted. qy: numpy.2darray Contains the y-coordinate values to be plotted qu: numpy.2darray Contains horizontal motion field values. qv: numpy.2darray Contains vertical motion field values. basetime: pandas.Timestamp Basetime of the nowcast. interval: pandas.Timedelta The interval between start and end time. coastline: cartopy.feature Feature (coastlines) to be added. """ cmap = plt.get_cmap(name='Reds') images = [] for t in ltgvf.coords['time'].values: plt.figure(figsize=(28, 21)) ax = plt.axes(projection=crs) ax.set_extent((712000, 962000, 695000, 945000), crs=crs) ax.add_feature(hires) ltgvf_t = ltgvf.sel(time=t) quadmesh = ltgvf.sel(time=t).plot(cmap=cmap) cbar = quadmesh.colorbar cbar.ax.set_ylabel( ltgvf.attrs['long_name'], fontsize=28 ) cbar.ax.tick_params(labelsize=24) ax.quiver( qx, qy, qu, qv, pivot='mid' ) ax.xaxis.set_visible(True) ax.yaxis.set_visible(True) ax.xaxis.set_tick_params(labelsize=24) ax.yaxis.set_tick_params(labelsize=24) ax.xaxis.label.set_size(28) ax.yaxis.label.set_size(28) basetime_str = basetime.strftime('%Y%m%d%H%M') validtime_str = pd.Timestamp(t).strftime('%Y%m%d%H%M') validtime = pd.Timestamp(t) t_minus = validtime-basetime-interval t_minus = int(t_minus.total_seconds()//60) t_plus = validtime-basetime t_plus = int(t_plus.total_seconds()//60) plt.title( f'Lightning Potential Forecast \nBasetime:{str(basetime)}' f' Valid time:{str(validtime)}\n' f'Start time: T{t_minus:+} minutes' f' End time: T{t_plus:+} minutes', fontsize=32 ) plt.savefig( THIS_DIR + "/../tests/outputs/ltgv_nowcast_" f"{basetime_str}_{validtime_str}.png" ) images.append( imageio.imread( THIS_DIR + "/../tests/outputs/ltgv_nowcast_" f"{basetime_str}_{validtime_str}.png" ) ) imageio.mimsave( THIS_DIR + f'/../tests/outputs/ltgv_nowcast.gif', images, duration=1 ) .. GENERATED FROM PYTHON SOURCE LINES 460-461 Step 3: Call the functions above to plot. .. GENERATED FROM PYTHON SOURCE LINES 461-511 .. code-block:: default # Coastlines hires = cfeature.GSHHSFeature( scale='h', levels=[1], edgecolor='black', facecolor='none' ) # Motion quivers # 256km radius qxres = area_def.x_size // 20 qyres = area_def.y_size // area_def.x_size * qxres qx = u.coords['easting'].values[::qxres] qy = u.coords['northing'].values[::qyres] qu = u.values[::qyres, ::qxres] qv = v.values[::qyres, ::qxres] # Zoomed in to Hong Kong uhk = u.sel( easting=slice(712000, 962000), northing=slice(945000, 695000) ) vhk = v.sel( easting=slice(712000, 962000), northing=slice(945000, 695000) ) qxreshk = qxres // 3*3 qyreshk = qxreshk // 2 quhk = uhk.values[::qyreshk, ::qxreshk] qvhk = vhk.values[::qyreshk, ::qxreshk] qxhk = uhk.coords['easting'].values[::qxreshk] qyhk = uhk.coords['northing'].values[::qyreshk] # projection crs = area_def.to_cartopy_crs() # Calling functions plot_reflectivity( reflectivity, crs, qx, qy, qu, qv, hires ) plot_nowcast( ltgvf, crs, qxhk, qyhk, quhk, qvhk, basetime, interval, hires ) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.003 seconds) .. _sphx_glr_download_auto_examples_ltg_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_hk.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ltg_hk.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_