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 time
from copy import copy
import pyproj
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm, ListedColormap, LinearSegmentedColormap
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_raw
from swirlspy.obs.lightning import Lightning
from swirlspy.ltg.map import gen_ltgv
from swirlspy.utils import standardize_attr, FrameType
plt.switch_backend('agg')
THIS_DIR = os.getcwd()
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.
# 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
)
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_geodetic = []
for idx, lst in enumerate(located_files):
ltg_object_geodetic = Lightning(
lst, 'geodetic', endtimes[idx]-interval, endtimes[idx]
)
ltg_object_list_geodetic.append(ltg_object_geodetic)
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_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.
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 = xr.concat(reflec_list, dim='time')
standardize_attr(reflectivity, frame_type=FrameType.dBZ)
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_geodetic,
area_def
):
"""
Convert coordinates in object from geodetic to HK1980.
Parameters
-----------
ltg_object_geodetic: swirlspy.obs.lightning.Lightning
Lightning object with geodetic 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_geodetic.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_geodetic)
ltg_object.lx = lx
ltg_object.ly = ly
ltg_object.data = data
ltg_object.proj = 'HK1980'
return ltg_object
# Calling the function to convert projection system of coordinates
ltg_object_list = []
for ltg_object_geodetic in ltg_object_list_geodetic:
ltg_object = convert_ltg(ltg_object_geodetic, 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 = xr.concat(ltgv_list, dim='time')
standardize_attr(ltgv, frame_type=FrameType.dBZ)
Generating motion fields and extrapolating¶
- Generate motion field using ROVER.
- Extrapolate lightning potential field using the motion field by Semi-Lagrangian Advection
# ROVER
motion = rover(reflectivity)
# SLA
steps = int(nowcast_duration.total_seconds()) // int(interval.total_seconds())
ltgvf = sla(ltgv, motion, nowcast_steps=steps)
Out:
RUNNING 'rover' FOR EXTRAPOLATION.....
Plotting results¶
Plot motion fields with reflectivities and lightning potential nowcasts.
Plotting reflectivities with motion fields.
# 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 = 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 = motion.coords['easting'].values[::5]
qy = motion.coords['northing'].values[::5]
qu = motion.values[0, ::5, ::5]
qv = motion.values[1, ::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
)
Plotting lightning potentials. In this example, only hourly potentials are plotted.
# 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
motion_sel = motion.sel(
easting=slice(zoom[0], zoom[1]),
northing=slice(zoom[3], zoom[2])
)
qx = motion_sel.coords['easting'].values[::10]
qy = motion_sel.coords['northing'].values[::10]
qu = motion_sel.values[0, ::10, ::10]
qv = motion_sel.values[1, ::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
)
Total running time of the script: ( 0 minutes 55.995 seconds)