{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Blending\nThis example demonstrates how to blend different \nreflectivity sources into one.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definitions\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\nimport matplotlib.pyplot as plt\nimport cartopy.crs as ccrs\nimport cartopy.feature as cfeature\nimport cartopy.io.shapereader as shpreader\nfrom matplotlib.gridspec import GridSpec\nfrom matplotlib.colors import BoundaryNorm\nfrom matplotlib.colors import ListedColormap\nfrom matplotlib.cm import ScalarMappable\nfrom pyresample import utils\n\nfrom swirlspy.core.resample import grid_resample\nfrom swirlspy.rad.iris import read_iris_grid\nfrom swirlspy.sat.h8 import read_h8_data\nfrom swirlspy.blending import comp_qpe, Raw\n\nplt.switch_backend('agg')\n\nroot_dir = os.getcwd()\n\nstart_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialising\n\nThis section demonstrates parsing\nHimawari-8 data.\n\nStep 1: Define necessary parameter.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Define base time\nbase_time = pd.Timestamp(\"2019-07-31T07:00\")\n\n# Define data boundary in WGS84 (latitude)\nlatitude_from = 30.\nlatitude_to = 16.\nlongitude_from = 105.\nlongitude_to = 122.\n\narea = (\n latitude_from, latitude_to,\n longitude_from, longitude_to\n)\n\n# Define grid size, use negative value for descending range\ngrid_size = (-.025, .025)\n\n# list of source data\nsources = []\n\ninitialising_time = pd.Timestamp.now()\n\n# Load map shape \nmap_shape_file = os.path.abspath(os.path.join(\n root_dir,\n '../tests/samples/shape/se_asia'\n))\n\n# coastline and province\nmap_with_province = cfeature.ShapelyFeature(\n list(shpreader.Reader(map_shape_file).geometries()),\n ccrs.PlateCarree()\n)\n\n# define the plot function\ndef plot_base(ax: plt.Axes, extents: list, crs: ccrs.Projection):\n ax.set_extent(extents, crs=crs)\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=ccrs.PlateCarree(),\n extent=[-180, 180, -180, 180], 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": [ "Step 2: Read data from radar files into xarray.DataArray\nusing read_iris_grid().\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "radar = read_iris_grid(\n os.path.join(root_dir, \"../tests/samples/iris/RAD190731150000.REF2256\")\n)\n\nradar_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Define the target grid as a pyresample AreaDefinition.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining target grid\narea_id = \"WGS84\"\ndescription = 'World Geodetic System 1984'\nproj_id = 'WGS84'\nprojection = '+proj=longlat +datum=WGS84 +no_defs'\nx_size = (longitude_to - longitude_from) / grid_size[1] + 1\ny_size = (latitude_to - latitude_from) / grid_size[0] + 1\narea_extent = (longitude_from, latitude_from, longitude_to, latitude_to)\nradar_area_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 5: Reproject the radar data from read_iris_grid() from Centered\nAzimuthal (source) projection to World Geodetic System 1984 projection.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Extracting the AreaDefinition of the source projection\narea_def_src = radar.attrs['area_def']\n\n# Reprojecting\nreproj_radar = grid_resample(\n radar, area_def_src, radar_area_def,\n coord_label=['x', 'y']\n).sortby(\n ['y'], ascending=False\n)\n\n# fix floating point issue\ny_coords = np.linspace(\n latitude_from,\n latitude_to,\n reproj_radar.data.shape[1],\n dtype=np.float32\n)\nx_coords = np.linspace(\n longitude_from,\n longitude_to,\n reproj_radar.data.shape[2],\n dtype=np.float32\n)\nreproj_radar.coords['y'] = np.array(y_coords)\nreproj_radar.coords['x'] = np.array(x_coords)\nreproj_radar = reproj_radar.sel(time=reproj_radar.coords['time'].values[0])\n\n\nradar_site = (\n reproj_radar.attrs['proj_site'][1],\n reproj_radar.attrs['proj_site'][0],\n 1.8, # radius\n 0.76 # weight sigma\n)\n\nsources.append(Raw(\n reproj_radar,\n [radar_site], # sites configuration, list of available sites useful for mosaic data\n 0.1 # data weight\n))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 6: Define data directory\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Supply data directory.\n# Please make sure H8 data filename is follow the naming pattern -\n# HS_H08_{date}_{time}_B{channel:02}_FLDK_R{rsol:02}_S{seg:02}10.DAT\n# example:\n# base time = 2019-07-31 07:00 UTC\n# channel = 4\n# resolution = 10\n# segment = 2\n# ========================================\n# filename: HS_H08_20190731_0700_B04_FLDK_R10_S0410.DAT\ndata_dir = os.path.join(root_dir, \"../tests/samples/h8\")\n\nsat_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 7: Parse data into reflectivity as xarray.DataArray\nusing read_h8_data().\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sat = read_h8_data(\n data_dir,\n base_time,\n area,\n grid_size\n)\n# remove time axis\nsat = sat.sel(time=sat.coords['time'].values[0])\n\n# no site data used, treat all points of data with same weight\nsources.append(Raw(\n sat,\n weight=0.01 # data weight\n))\n\nblend_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 8: Blend all data together.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflec = comp_qpe(\n grid_size,\n area,\n sources\n)\n\n\npost_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 9: Remove invalid data if needed.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflec.values[reflec.values < 13.] = reflec.attrs['zero_value']\n\n# update sat data for plotting\nsat.values[sat.values < 13.] = sat.attrs['zero_value']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating radar reflectivity maps\n\nDefine the color scale and format of the plots\nand plot using xarray.plot().\n\nIn this example, only hourly images will be plotted.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining colour scale and format\nlevels = [\n -32768,\n 10, 15, 20, 24, 28, 32,\n 34, 38, 41, 44, 47, 50,\n 53, 56, 58, 60, 62\n]\ncmap = ListedColormap([\n '#FFFFFF00', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433',\n '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200',\n '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0'\n])\n\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\n# colorbar map\nmappable = ScalarMappable(cmap=cmap, norm=norm)\nmappable.set_array([])\n\n# Defining the crs\ncrs = ccrs.PlateCarree()\nextents = (longitude_from, longitude_to, latitude_from, latitude_to)\n\n# Plotting\nfig: plt.Figure = plt.figure(figsize=(24, 8), frameon=False)\ngs = GridSpec(\n 1, 3, figure=fig, wspace=0.03, hspace=-0.25, top=0.95,\n bottom=0.05, left=0.17, right=0.845\n)\n\n# plot radar\nax = fig.add_subplot(gs[0, 0], projection=crs)\n\nplot_base(ax, extents, crs)\n\nim = ax.imshow(reproj_radar.values, cmap=cmap, norm=norm, interpolation='nearest',\n extent=extents)\n\nax.set_title(\n \"RADAR\\n\"\n f\"Based @ {base_time.strftime('%H:%MH')}\",\n loc='left',\n fontsize=9\n)\nax.set_title(\n ''\n)\nax.set_title(\n f\"{base_time.strftime('%Y-%m-%d')} \\n\"\n f\"Valid @ {(base_time + pd.Timedelta(minutes=10)).strftime('%H:%MH')} \",\n loc='right',\n fontsize=9\n)\n\n# plot H8\nax = fig.add_subplot(gs[0, 1], projection=crs)\n\nplot_base(ax, extents, crs)\n\nim = ax.imshow(sat.values, cmap=cmap, norm=norm, interpolation='nearest',\n extent=extents)\n\nax.set_title(\n \"H8\\n\"\n f\"Based @ {base_time.strftime('%H:%MH')}\",\n loc='left',\n fontsize=9\n)\nax.set_title(\n ''\n)\nax.set_title(\n f\"{base_time.strftime('%Y-%m-%d')} \\n\"\n f\"Valid @ {(base_time + pd.Timedelta(minutes=10)).strftime('%H:%MH')} \",\n loc='right',\n fontsize=9\n)\n\n# plot blended\nax = fig.add_subplot(gs[0, 2], projection=crs)\n\nplot_base(ax, extents, crs)\n\nim = ax.imshow(reflec.values, cmap=cmap, norm=norm, interpolation='nearest',\n extent=extents)\n\nax.set_title(\n \"Blended\\n\"\n f\"Based @ {base_time.strftime('%H:%MH')}\",\n loc='left',\n fontsize=9\n)\nax.set_title(\n ''\n)\nax.set_title(\n f\"{base_time.strftime('%Y-%m-%d')} \\n\"\n f\"Valid @ {(base_time + pd.Timedelta(minutes=10)).strftime('%H:%MH')} \",\n loc='right',\n fontsize=9\n)\n\ncbar_ax = fig.add_axes([0.875, 0.125, 0.03, 0.75])\ncbar = fig.colorbar(\n mappable, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g')\ncbar.ax.set_ylabel('Reflectivity', rotation=90)\n\nfig.savefig(\n os.path.join(\n root_dir,\n \"../tests/outputs/blending.png\"\n ),\n bbox_inches='tight'\n)\n\nimage_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Checking run time of each component\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(f\"Start time: {start_time}\")\nprint(f\"Initialising time: {initialising_time}\")\nprint(f\"Read radar time: {radar_time}\")\nprint(f\"Parse H8 data time: {sat_time}\")\nprint(f\"Blending time: {blend_time}\")\nprint(f\"Post blending time: {post_time}\")\nprint(f\"Plotting blended image time: {image_time}\")\n\nprint(f\"Time to initialise: {initialising_time - start_time}\")\nprint(f\"Time to run read radar: {radar_time - initialising_time}\")\nprint(f\"Time to run data parsing: {sat_time - radar_time}\")\nprint(f\"Time to run blending: {blend_time - sat_time}\")\nprint(f\"Time to perform post process: {post_time - blend_time}\")\nprint(f\"Time to plot reflectivity image: {image_time - post_time}\")" ] } ], "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 }