Note
Click here to download the full example code
Lightning Nowcast (Hong Kong)
This example demonstrates how to perform lightning nowcasts of up to three hours, using data from Hong Kong.
Definitions
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()
Traceback (most recent call last):
File "/tmp/build/docs/swirlspy/swirlspy/examples/ltg_hk.py", line 28, in <module>
from swirlspy.qpf.rover import rover
File "/tmp/build/docs/swirlspy/swirlspy/qpf/rover.py", line 6, in <module>
from rover.rover import rover as rover_api
ImportError: libopencv_core.so.3.4: cannot open shared object file: No such file or directory
Initialising
Step 1: Define the target grid as a pyresample AreaDefinition.
# 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
)
Step 2: Defining a basetime, nowcast interval and nowcast duration.
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)
Step 3: Extracting Lightning Location Information System (LLIS) lightning files to read.
# 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)
Step 4: Read data from files into Lightning objects.
# 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)
Step 5: Locating IRIS RAW radar files needed to generate motion field.
# 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))
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.
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')
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.
# 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
# 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)
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.
# 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')
Generating motion fields and extrapolating
Generate motion field using ROVER.
Extrapolate lightning potential field using the motion field by Semi-Lagrangian Advection
# ROVER
u, v = rover(reflectivity)
# SLA
steps = int(nowcast_duration.total_seconds()) // int(interval.total_seconds())
ltgvf = sla(ltgv, u, v, steps=steps)
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.
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
)
Step 2: Write a function to plot nowcasted lightning potential.
# 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
)
Step 3: Call the functions above to plot.
# 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
)
Total running time of the script: ( 0 minutes 0.003 seconds)