{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nSimple Lightning Nowcast (Hong Kong)\n========================================================\nThis example demonstrates how to perform\na simple lightning nowcast by the extrapolation\nof lightning strikes, using data from Hong Kong.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Definitions\n--------------------------------------------------------\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import all required modules and methods:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Python package to allow system command line functions\nimport os\n# Python package to manage warning message\nimport warnings\n# Python package for time calculations\nimport pandas as pd\n# Python package for numerical calculations\nimport numpy as np\n# Python package for xarrays to read and handle netcdf data\nimport xarray as xr\n# Python package for text formatting\nimport textwrap\n# Python package for projection description\nimport pyproj\nfrom pyresample import get_area_def\n# Python package for projection\nimport cartopy.crs as ccrs\n# Python package for land/sea features\nimport cartopy.feature as cfeature\n# Python package for reading map shape file\nimport cartopy.io.shapereader as shpreader\n# Python package for creating plots\nfrom matplotlib import pyplot as plt\n# Python package for output grid format\nfrom matplotlib.gridspec import GridSpec\n# Python package for colorbars\nfrom matplotlib.colors import BoundaryNorm, ListedColormap\nfrom matplotlib.cm import ScalarMappable\n\n# Python com-swirls package to calculate motion field (rover) and semi-lagrangian advection\nfrom swirlspy.qpf import rover, sla\n# swirlspy data parser function\nfrom swirlspy.rad.iris import read_iris_grid\n# swirlspy test data source locat utils function\nfrom swirlspy.qpe.utils import timestamps_ending, locate_file\n# swirlspy regrid function\nfrom swirlspy.core.resample import grid_resample\n# swirlspy lightning function\nfrom swirlspy.obs import Lightning\nfrom swirlspy.ltg.map import gen_ltg_binary\n# swirlspy standardize data function\nfrom swirlspy.utils import standardize_attr, FrameType\n# directory constants\nfrom swirlspy.tests.samples import DATA_DIR\nfrom swirlspy.tests.outputs import OUTPUT_DIR\n\nwarnings.filterwarnings(\"ignore\")\nplt.switch_backend('agg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initialising\n-------------------------------------------------------------------------\n\nDefine the target grid:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining target grid\narea_id = \"hk1980_250km\"\ndescription = (\"A 1km resolution rectangular grid \"\n \"centred at HKO and extending to 250 km \"\n \"in each direction in HK1980 easting/northing coordinates\")\nproj_id = 'hk1980'\nprojection = ('+proj=tmerc +lat_0=22.31213333333334 '\n '+lon_0=114.1785555555556 +k=1 +x_0=836694.05 '\n '+y_0=819069.8 +ellps=intl +towgs84=-162.619,-276.959,'\n '-161.764,0.067753,-2.24365,-1.15883,-1.09425 +units=m '\n '+no_defs')\n\nx_size = 500\ny_size = 500\n\narea_extent = (587000, 569000, 1087000, 1069000)\n\narea_def = get_area_def(\n area_id, description, proj_id, projection, x_size, y_size, area_extent\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define basemap:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Load the shape of Hong Kong\nmap_shape_file = os.path.abspath(os.path.join(\n DATA_DIR,\n 'shape/hk_hk1980'\n))\n\n# coastline and province\nmap_with_province = cfeature.ShapelyFeature(\n list(shpreader.Reader(map_shape_file).geometries()),\n area_def.to_cartopy_crs()\n)\n\n\n# define the plot function\ndef plot_base(ax: plt.Axes, extents: list, crs: ccrs.Projection):\n # fake the ocean color\n ax.imshow(\n np.tile(np.array([[[178, 208, 254]]], dtype=np.uint8), [2, 2, 1]),\n origin='upper', transform=crs,\n extent=extents, zorder=-1\n )\n # coastline, province and state, color\n ax.add_feature(\n map_with_province, facecolor=cfeature.COLORS['land'],\n edgecolor='none', zorder=0\n )\n # overlay coastline, province and state without color\n ax.add_feature(\n map_with_province, facecolor='none', edgecolor='gray', linewidth=0.5\n )\n ax.set_title('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Defining a basetime, nowcast interval and nowcast duration:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "basetime = pd.Timestamp('201904201200').floor('min')\nbasetime_str = basetime.strftime('%Y%m%d%H%M')\nnowcast_interval = pd.Timedelta(30, 'm')\nnowcast_duration = pd.Timedelta(2, 'h')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Loading Lightning and Radar Data\n--------------------------------------------------------------------------\n\nExtract Lightning Location Information System (LLIS) files to read.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# The interval between LLIS files\nllis_file_interval = pd.Timedelta(1, 'm')\n\n# Finding the string representation\n# of timestamps marking the relevant LLIS files\ntimestamps = timestamps_ending(\n basetime,\n duration=nowcast_interval,\n interval=llis_file_interval,\n format='%Y%m%d%H%M'\n)\n\n# Locating files\nllis_dir = os.path.join(DATA_DIR, 'llis')\n\nlocated_files = []\nfor timestamp in timestamps:\n located_file = locate_file(llis_dir, timestamp)\n located_files.append(located_file)\n if located_file is not None:\n print(f\"LLIS file for {timestamp} has been located\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read data from files into Lightning objects.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Initialising Lightning objects\n# Coordinates are geodetic\nltg_object_geodetic = Lightning(\n located_files,\n 'geodetic',\n basetime - nowcast_interval,\n basetime\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Locating IRIS REF radar files needed to generate motion field.\nThe IRIS REF files are marked by start time.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# The interval between the radar files\nrad_interval = pd.Timedelta(6, 'm')\n\n# Finding the string representation\n# of the timestamps\n# marking the required files\n# In this case, the most recent radar files\n# separated by `nowcast_interval` is required\n\nrad_timestamps = timestamps_ending(\n basetime,\n duration=nowcast_interval + rad_interval,\n interval=rad_interval\n)\n\nrad_timestamps = [rad_timestamps[0], rad_timestamps[-1]]\n\nrad_dir = os.path.join(DATA_DIR, 'iris')\nlocated_rad_files = []\nfor timestamp in rad_timestamps:\n located_file = locate_file(rad_dir, timestamp)\n located_rad_files.append(located_file)\n if located_file is not None:\n print(f\"IRIS REF file for {timestamp} has been located\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read from iris radar files using read_iris_grid().\nThe data is in the Azimuthal Equidistant Projection.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflec_unproj_list = []\nfor file in located_rad_files:\n reflec = read_iris_grid(file)\n reflec_unproj_list.append(reflec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Data Reprojection\n---------------------------------------------------------------------------\n\nSince the coordinates of the Lightning objects are geodetic, a conversion\nfrom geodetic coordinates to that of the user-defined projection system is\nrequired.\nThis can be achieved by `swirlspy.obs.Lightning.reproject()`\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# pyproj representation of the map projection of the input data\ninProj = pyproj.Proj(init='epsg:4326')\n# pyproj representation of the intended map projection\noutProj = pyproj.Proj(area_def.proj_str)\n\n# reproject\nltg_object = ltg_object_geodetic.reproject(inProj, outProj, 'HK1980')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Radar data also needs to be reprojected from Azimuthal Equidistant\nto HK1980. This is achieved by the `swirlspy.core.resample` function.\nFollowing reprojection, the individual xarray.DataArrays can be concatenated\nto a single xarray.DataArrays.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflec_list = []\nfor z in reflec_unproj_list:\n reflec = grid_resample(\n z,\n z.attrs['area_def'],\n area_def,\n coord_label=['easting', 'northing']\n )\n reflec_list.append(reflec)\n\nreflectivity = xr.concat(reflec_list, dim='time')\nstandardize_attr(reflectivity, frame_type=FrameType.dBZ)\n\nprint(reflectivity)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Grid lightning observations\n-----------------------------------------------------------------\n\nGenerate a binary grid of lighting observations, where 0\nindicates no lightning and 1 indicates that there is lightning.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ltg_bin = gen_ltg_binary(ltg_object,\n area_def,\n coord_label=['easting', 'northing'])\nprint(ltg_bin)\n\n# Duplicating ltg_density along the time dimension\n# so that it can be advected.\nltg_bin1 = ltg_bin.copy()\nltg_bin1.coords['time'] = [basetime - nowcast_interval]\nltg_bin_concat = xr.concat(\n [ltg_bin1, ltg_bin],\n dim='time'\n)\n# Identifying zero value\nltg_bin_concat.attrs['zero_value'] = np.nan" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generating motion fields and extrapolating\n--------------------------------------------------------------\n\n#. Generate motion field using ROVER.\n#. Extrapolate gridded lightning strikes using the motion field\n by Semi-Lagrangian Advection.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "motion = rover(reflectivity) # rover\n\nsteps = int(nowcast_duration.total_seconds()\n // nowcast_interval.total_seconds())\nltg_frames = sla(ltg_bin_concat, motion, nowcast_steps=steps-1) # sla\nprint(ltg_frames)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualization\n---------------------------------------------------------------\n\nNow, let us visualize the lightning nowcast up to a duration\nof 2 hours.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Colormap and colorbar parameters\nn_col = 2\ncmap = ListedColormap(['#FFFFFF', '#000000'])\nlevels = [0., 0.5, 1.]\nticks = np.array([0, 1])\nnorm = BoundaryNorm(levels, cmap.N, clip=True)\n\nmappable = ScalarMappable(cmap=cmap, norm=norm)\nmappable.set_array([])\n\n# Defining the crs\ncrs = area_def.to_cartopy_crs()\n\n# Defining area\nx = ltg_frames.coords['easting'].values\ny = ltg_frames.coords['northing'].values\nextents = [\n np.min(x), np.max(x),\n np.min(y), np.max(y)\n]\n\nqx = motion.coords['easting'].values[::20]\nqy = motion.coords['northing'].values[::20]\nqu = motion.values[0, ::20, ::20]\nqv = motion.values[1, ::20, ::20]\n\nfig = plt.figure(figsize=(8, 8), frameon=False)\ngs = GridSpec(\n 2, 2, figure=fig, wspace=0.03, hspace=-0.25, top=0.95,\n bottom=0.05, left=0.17, right=0.845\n)\n\nfor i, t in enumerate(ltg_frames.time.values):\n time = pd.Timestamp(t)\n\n row = i // 2\n col = i % 2\n ax = fig.add_subplot(gs[row, col], projection=crs)\n\n # plot base map\n plot_base(ax, extents, crs)\n\n # plot lightning\n frame = ltg_frames.sel(time=time)\n frame.where(frame > levels[1]).plot(\n cmap=cmap,\n norm=norm,\n add_colorbar=False\n )\n\n # plot motion vector\n ax.quiver(\n qx, qy, qu, qv, pivot='mid'\n )\n\n ax.set_title('')\n\n ax.text(\n extents[0],\n extents[3],\n textwrap.dedent(\n \"\"\"\n Lightning Density\n Based @ {baseTime}\n \"\"\"\n ).format(\n baseTime=basetime.strftime('%H:%MH')\n ).strip(),\n fontsize=10,\n va='bottom',\n ha='left',\n linespacing=1\n )\n\n ax.text(\n extents[1],\n extents[3],\n textwrap.dedent(\n \"\"\"\n {validDate}\n Valid @ {validTime}\n \"\"\"\n ).format(\n validDate=basetime.strftime('%Y-%m-%d'),\n validTime=time.strftime('%H:%MH')\n ).strip(),\n fontsize=10,\n va='bottom',\n ha='right',\n linespacing=1\n )\n\ncbar_ax = fig.add_axes([0.875, 0.125, 0.03, 0.75])\ncbar = fig.colorbar(mappable, cax=cbar_ax)\ntick_locs = (np.arange(n_col) + 0.5)*(n_col-1)/n_col\ncbar.set_ticks(tick_locs)\ncbar.set_ticklabels(np.arange(n_col))\ncbar.ax.set_ylabel(\n f\"{ltg_bin.attrs['long_name']}[{ltg_bin.attrs['units']}]\"\n)\n\nfig.savefig(\n os.path.join(OUTPUT_DIR, \"ltg-bin-hk.png\"),\n bbox_inches='tight'\n)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.15" } }, "nbformat": 4, "nbformat_minor": 0 }