.. 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_raw from swirlspy.obs.lightning import Lightning from swirlspy.ltg.map import gen_ltgv from swirlspy.core.resample import grid_resample plt.switch_backend('agg') 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-47 Initialising ------------------------------------------------------------------------- Step 1: Define the target grid as a pyresample AreaDefinition. If a zoom in during the plotting phase is required, a separate AreaDefinition can be defined. .. GENERATED FROM PYTHON SOURCE LINES 47-75 .. 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 ) area_extent_plot = (712000, 695000, 962000, 945000) area_def_plot = utils.get_area_def( area_id, description, proj_id, projection, x_size, y_size, area_extent_plot ) .. GENERATED FROM PYTHON SOURCE LINES 76-78 Step 2: Defining a basetime, nowcast interval and nowcast duration. .. GENERATED FROM PYTHON SOURCE LINES 78-85 .. 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 86-89 Step 3: Extracting Lightning Location Information System (LLIS) lightning files to read. .. GENERATED FROM PYTHON SOURCE LINES 89-112 .. 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 113-115 Step 4: Read data from files into Lightning objects. .. GENERATED FROM PYTHON SOURCE LINES 115-125 .. 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 126-127 Step 5: Locating IRIS RAW radar files needed to generate motion field. .. GENERATED FROM PYTHON SOURCE LINES 127-138 .. 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 139-143 Step 6: Read from iris radar files using read_iris_raw(). Projection system is specified in the read_iris_raw() function using the AreaDefinition defined above. The resulting xarray.DataArrays are then concatenated along the time axis. .. GENERATED FROM PYTHON SOURCE LINES 143-154 .. code-block:: default reflec_list = [] for file in located_rad_files: reflec = read_iris_raw( 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 155-161 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 161-213 .. 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 214-223 .. 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 224-232 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 232-250 .. 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 251-258 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 258-266 .. 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 267-272 Plotting results ---------------------------------------------------------------- Plot motion fields with reflectivities and lightning potential nowcasts. .. GENERATED FROM PYTHON SOURCE LINES 274-275 Plotting reflectivities with motion fields. .. GENERATED FROM PYTHON SOURCE LINES 275-337 .. code-block:: default # Defining colour scale 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) # Defining the crs crs = area_def.to_cartopy_crs() qx = u.coords['easting'].values[::5] qy = u.coords['northing'].values[::5] qu = u.values[::5, ::5] qv = v.values[::5, ::5] # Plotting p = reflectivity.plot( col='time', subplot_kws={'projection': crs}, cbar_kwargs={ 'extend': 'max', 'ticks': levels[1:], 'format': '%.3g' }, cmap=cmap, norm=norm ) for idx, ax in enumerate(p.axes.flat): time = pd.Timestamp(reflectivity.coords['time'].values[idx]) ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20) ax.coastlines('10m') ax.gridlines() ax.set_title( "Reflectivity\n" f"Based @ {basetime.strftime('%H:%MH')}", loc='left', fontsize=9 ) ax.set_title( '' ) ax.set_title( f"{basetime.strftime('%Y-%m-%d')} \n" f"Valid @ {time.strftime('%H:%MH')} ", loc='right', fontsize=9 ) plt.savefig( THIS_DIR + f"/../tests/outputs/reflectivity-hk.png", dpi=300 ) .. GENERATED FROM PYTHON SOURCE LINES 338-341 Plotting lightning potentials. In this example, only hourly potentials are plotted. .. GENERATED FROM PYTHON SOURCE LINES 341-413 .. code-block:: default # Defining the crs crs = area_def_plot.to_cartopy_crs() # Defining zoom zoom = (712000, 962000, 695000, 945000) # (x0, x1, y0, y1) # Generating a timelist for every hour timelist = [ basetime + pd.Timedelta(minutes=60*i) for i in range(4) ] # Obtaining the slice of the xarray to be plotted da_plot = ltgvf.sel( time=timelist, easting=slice(zoom[0], zoom[1]), northing=slice(zoom[3], zoom[2]) ) cmap = plt.get_cmap('Reds') # Defining motion quivers u_sel = u.sel( easting=slice(zoom[0], zoom[1]), northing=slice(zoom[3], zoom[2]) ) v_sel = v.sel( easting=slice(zoom[0], zoom[1]), northing=slice(zoom[3], zoom[2]) ) qx = u_sel.coords['easting'].values[::10] qy = u_sel.coords['northing'].values[::10] qu = u_sel.values[::10, ::10] qv = v_sel.values[::10, ::10] # Plotting p = da_plot.plot( col='time', col_wrap=2, subplot_kws={'projection': crs}, cbar_kwargs={ 'extend': 'max', 'ticks': [i/10 for i in range(11)], 'format': '%.3g' }, cmap=cmap ) for idx, ax in enumerate(p.axes.flat): ax.quiver(qx, qy, qu, qv, pivot='mid') ax.coastlines('10m') ax.set_xlim(zoom[0], zoom[1]) ax.set_ylim(zoom[2], zoom[3]) ax.gridlines() ax.set_title( "Lightning Potential\n" f"Based @ {basetime.strftime('%H:%MH')}", loc='left', fontsize=8 ) ax.set_title( '' ) ax.set_title( f"{basetime.strftime('%Y-%m-%d')} \n" f"Valid @ {timelist[idx].strftime('%H:%MH')} ", loc='right', fontsize=8 ) plt.savefig( THIS_DIR + f"/../tests/outputs/ltgvf.png", dpi=300 ) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.017 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 `_