{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Lightning Nowcast (Hong Kong)\nThis example demonstrates how to perform\nlightning nowcasts of up to three hours, using data\nfrom Hong Kong.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definitions\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nimport xarray\nimport time\nimport pandas as pd\nimport matplotlib\nimport imageio\nimport pyproj\nimport matplotlib.pyplot as plt\nimport cartopy.feature as cfeature\nimport cartopy.crs as ccrs\nfrom pyresample import utils\nfrom copy import copy\nfrom matplotlib.colors import LinearSegmentedColormap, BoundaryNorm\n\nfrom swirlspy.qpf.rover import rover\nfrom swirlspy.qpf.sla import sla\nfrom swirlspy.qpe.utils import timestamps_ending, locate_file\nfrom swirlspy.rad.iris import read_iris\nfrom swirlspy.obs.lightning import Lightning\nfrom swirlspy.ltg.map import gen_ltgv\nfrom swirlspy.core.resample import grid_resample\n\n\nTHIS_DIR = os.getcwd()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialising\n\nStep 1: Define the target grid as a pyresample AreaDefinition.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining target grid\narea_id = \"hk1980_250km\"\ndescription = (\"A 250 m 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 = utils.get_area_def(\n area_id, description, proj_id, projection, x_size, y_size, area_extent\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: 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('20180811073000').floor('min')\nbasetime_str = basetime.strftime('%Y%m%d%H%M')\ninterval = pd.Timedelta(minutes=15)\nnowcast_duration = pd.Timedelta(hours=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Extracting Lightning Location Information System (LLIS)\nlightning files to read.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining the times of files and the accumulation endtimes of files\nendtimes = [basetime - interval, basetime]\ntimestamps = []\nfor time in endtimes:\n ts = timestamps_ending(\n time,\n duration=interval,\n interval=pd.Timedelta(minutes=1)\n )\n timestamps.append(ts)\n\n# Locating files\nllis_dir = THIS_DIR + '/../tests/samples/llis'\nlocated_files = []\n\nfor tss in timestamps:\n lf = []\n for timestamp in tss:\n lf.append(locate_file(llis_dir, timestamp))\n located_files.append(lf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 4: 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_list_wgs84 = []\nfor idx, lst in enumerate(located_files):\n ltg_object_wgs84 = Lightning(\n lst, 'WGS84', endtimes[idx]-interval, endtimes[idx]\n )\n ltg_object_list_wgs84.append(ltg_object_wgs84)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 5: Locating IRIS RAW radar files needed to generate motion field.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Generating timestrings for required times\nendtimes_str = []\nfor time in endtimes:\n endtimes_str.append(time.strftime('%y%m%d%H%M'))\n\nrad_dir = THIS_DIR + '/../tests/samples/iris/'\nlocated_rad_files = []\nfor time_str in endtimes_str:\n located_rad_files.append(locate_file(rad_dir, time_str))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 6: Read from iris radar files using read_iris(). Projection system\nis specified in the read_iris() function using the AreaDefinition\ndefined above. The resulting xarray.DataArrays are then concatenated\nalong the time axis.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflec_list = []\nfor file in located_rad_files:\n reflec = read_iris(\n file, area_def=area_def,\n coord_label=['easting', 'northing']\n )\n reflec_list.append(reflec)\n\nreflectivity = xarray.concat(reflec_list, dim='time')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data Reprojection\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.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Define a function for neater conversion\n\n\ndef convert_ltg(\n ltg_object_wgs84,\n area_def\n):\n\n \"\"\"\n Convert coordinates in object from WGS84 to HK1980.\n\n Parameters\n -----------\n ltg_object_wgs84: swirlspy.obs.lightning.Lightning\n Lightning object with WGS84 coordinates.\n area_def: pyresample.geometry.AreaDefinition\n AreaDefinition of the target grid.\n\n Returns\n -------\n ltg_object: swirlspy.obs.lightning.Lightning\n Rain object in HK1980.\n\n \"\"\"\n\n # Defining Proj objects\n geod = pyproj.Proj(init=\"epsg:4326\")\n outproj = pyproj.Proj(area_def.proj_str)\n\n lx = []\n ly = []\n data = []\n\n for tup in ltg_object_wgs84.data:\n lat = tup[1]\n lon = tup[2]\n easting, northing = pyproj.transform(\n geod, outproj, lon, lat\n )\n lx.append(easting)\n ly.append(northing)\n data.append((tup[0], northing, easting))\n\n ltg_object = copy(ltg_object_wgs84)\n ltg_object.lx = lx\n ltg_object.ly = ly\n ltg_object.data = data\n ltg_object.proj = 'HK80'\n\n return ltg_object" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Calling the function to convert projection system of coordinates\n\nltg_object_list = []\nfor ltg_object_wgs84 in ltg_object_list_wgs84:\n ltg_object = convert_ltg(ltg_object_wgs84, area_def)\n ltg_object_list.append(ltg_object)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating lightning potential field\n\nGenerate lightning potential fields using gen_ltgv(). Lightning\npotential decrease normally according to the provided sigma from\nthe site. Lightning potentials are then concatenated along the time\naxis.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Obtaining lightning potential field from lightning objects\nltgv_list = []\nfor ltg_object in ltg_object_list:\n ltgvi = gen_ltgv(\n ltg_object,\n area_def,\n ltg_object.start_time,\n ltg_object.end_time,\n sigma=20000,\n coord_label=['easting', 'northing']\n )\n ltgv_list.append(ltgvi)\n\n\n# Concatenating lightning potential along time dimension\nltgv = xarray.concat(ltgv_list, dim='time')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating motion fields and extrapolating\n\n1. Generate motion field using ROVER.\n2. Extrapolate lightning potential field using the motion field\n by Semi-Lagrangian Advection\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# ROVER\nu, v = rover(reflectivity)\n# SLA\nsteps = int(nowcast_duration.total_seconds()) // int(interval.total_seconds())\nltgvf = sla(ltgv, u, v, steps=steps)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting results\n\nPlot motion fields with reflectivities and lightning potential nowcasts.\nFunctions can be written to make plotting neater.\n\n\nStep 1: Write function to plot reflectivity with motion fields.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def plot_reflectivity(\n reflectivity, crs, qx, qy, qu, qv,\n coastline\n):\n \"\"\"\n Custom function for plotting reflectivities\n with motion fields.\n\n Parameters\n ------------------------------------------\n reflectivity: xarray.DataArray\n Contains data to be plotted\n crs: cartopy.crs.CRS\n Defines the projection system.\n qx: numpy.2darray\n Contains the x-coordinate values to\n be plotted.\n qy: numpy.2darray\n Contains the y-coordinate values to\n be plotted\n qu: numpy.2darray\n Contains horizontal motion field values.\n qv: numpy.2darray\n Contains vertical motion field values.\n coastline: cartopy.feature\n Feature (coastline) to be added.\n\n \"\"\"\n # Defining colour levels and format\n levels = [\n -32768,\n 10, 15, 20, 24, 28, 32,\n 34, 38, 41, 44, 47, 50,\n 53, 56, 58, 60, 62\n ]\n cmap = matplotlib.colors.ListedColormap([\n '#FFFFFF', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433',\n '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200',\n '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0'\n ])\n norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\n images = []\n for time in reflectivity.coords['time'].values:\n # Plotting\n plt.figure(figsize=(28, 21))\n ax = plt.axes(projection=crs)\n ax.add_feature(hires)\n quadmesh = reflectivity.sel(time=time).plot(\n cmap=cmap, norm=norm, extend='neither'\n )\n cbar = quadmesh.colorbar\n cbar.ax.set_ylabel(\n reflectivity.attrs['long_name']+'[' +\n reflectivity.attrs['units']+']',\n fontsize=28\n )\n cbar.ax.tick_params(labelsize=24)\n ax.quiver(\n qx, qy, qu, qv,\n pivot='mid')\n ax.xaxis.set_visible(True)\n ax.yaxis.set_visible(True)\n ax.xaxis.set_tick_params(labelsize=24)\n ax.yaxis.set_tick_params(labelsize=24)\n ax.xaxis.label.set_size(28)\n ax.yaxis.label.set_size(28)\n pdtime = pd.Timestamp(time)\n plt.title(\n \"Reflectivity Field with Motion Vectors\\n\"\n f\"Time: {pdtime}\",\n fontsize=25\n )\n plt.savefig(\n THIS_DIR +\n f\"/../tests/outputs/z_{pdtime.strftime('%Y%m%d%H%M')}.png\"\n )\n images.append(\n imageio.imread(\n THIS_DIR +\n f\"/../tests/outputs/z_{pdtime.strftime('%Y%m%d%H%M')}.png\"\n )\n )\n imageio.mimsave(\n THIS_DIR + \"/../tests/outputs/reflectivity.gif\",\n images,\n duration=1\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: Write a function to plot nowcasted lightning potential.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Plotting potential\n\n\ndef plot_nowcast(\n ltgvf, crs,\n qx, qy, qu, qv,\n basetime,\n interval,\n coastline\n):\n \"\"\"\n Custom function for plotting lightning potential nowcasts\n with motion fields.\n\n Parameters\n ------------------------------------------\n ltgvf: xarray.DataArray\n Contains data to be plotted\n crs: cartopy.crs.CRS\n Defines the projection system.\n qx: numpy.2darray\n Contains the x-coordinate values to\n be plotted.\n qy: numpy.2darray\n Contains the y-coordinate values to\n be plotted\n qu: numpy.2darray\n Contains horizontal motion field values.\n qv: numpy.2darray\n Contains vertical motion field values.\n basetime: pandas.Timestamp\n Basetime of the nowcast.\n interval: pandas.Timedelta\n The interval between start and end time.\n coastline: cartopy.feature\n Feature (coastlines) to be added.\n\n \"\"\"\n\n cmap = plt.get_cmap(name='Reds')\n images = []\n for t in ltgvf.coords['time'].values:\n plt.figure(figsize=(28, 21))\n ax = plt.axes(projection=crs)\n ax.set_extent((712000, 962000, 695000, 945000), crs=crs)\n ax.add_feature(hires)\n ltgvf_t = ltgvf.sel(time=t)\n quadmesh = ltgvf.sel(time=t).plot(cmap=cmap)\n cbar = quadmesh.colorbar\n cbar.ax.set_ylabel(\n ltgvf.attrs['long_name'],\n fontsize=28\n )\n cbar.ax.tick_params(labelsize=24)\n ax.quiver(\n qx, qy, qu, qv, pivot='mid'\n )\n ax.xaxis.set_visible(True)\n ax.yaxis.set_visible(True)\n ax.xaxis.set_tick_params(labelsize=24)\n ax.yaxis.set_tick_params(labelsize=24)\n ax.xaxis.label.set_size(28)\n ax.yaxis.label.set_size(28)\n\n basetime_str = basetime.strftime('%Y%m%d%H%M')\n validtime_str = pd.Timestamp(t).strftime('%Y%m%d%H%M')\n validtime = pd.Timestamp(t)\n t_minus = validtime-basetime-interval\n t_minus = int(t_minus.total_seconds()//60)\n t_plus = validtime-basetime\n t_plus = int(t_plus.total_seconds()//60)\n plt.title(\n f'Lightning Potential Forecast \\nBasetime:{str(basetime)}'\n f' Valid time:{str(validtime)}\\n'\n f'Start time: T{t_minus:+} minutes'\n f' End time: T{t_plus:+} minutes',\n fontsize=32\n )\n plt.savefig(\n THIS_DIR +\n \"/../tests/outputs/ltgv_nowcast_\"\n f\"{basetime_str}_{validtime_str}.png\"\n )\n\n images.append(\n imageio.imread(\n THIS_DIR +\n \"/../tests/outputs/ltgv_nowcast_\"\n f\"{basetime_str}_{validtime_str}.png\"\n )\n )\n imageio.mimsave(\n THIS_DIR +\n f'/../tests/outputs/ltgv_nowcast.gif',\n images,\n duration=1\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Call the functions above to plot.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Coastlines\nhires = cfeature.GSHHSFeature(\n scale='h',\n levels=[1],\n edgecolor='black',\n facecolor='none'\n )\n\n# Motion quivers\n# 256km radius\nqxres = area_def.x_size // 20\nqyres = area_def.y_size // area_def.x_size * qxres\n\nqx = u.coords['easting'].values[::qxres]\nqy = u.coords['northing'].values[::qyres]\nqu = u.values[::qyres, ::qxres]\nqv = v.values[::qyres, ::qxres]\n\n# Zoomed in to Hong Kong\nuhk = u.sel(\n easting=slice(712000, 962000), northing=slice(945000, 695000)\n )\nvhk = v.sel(\n easting=slice(712000, 962000), northing=slice(945000, 695000)\n )\nqxreshk = qxres // 3*3\nqyreshk = qxreshk // 2\n\nquhk = uhk.values[::qyreshk, ::qxreshk]\nqvhk = vhk.values[::qyreshk, ::qxreshk]\nqxhk = uhk.coords['easting'].values[::qxreshk]\nqyhk = uhk.coords['northing'].values[::qyreshk]\n\n# projection\ncrs = area_def.to_cartopy_crs()\n\n# Calling functions\nplot_reflectivity(\n reflectivity, crs,\n qx, qy, qu, qv,\n hires\n)\nplot_nowcast(\n ltgvf, crs,\n qxhk, qyhk, quhk, qvhk,\n basetime,\n interval,\n hires\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 }