{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# SAFNWC data\nThis example demonstrates how to read SAFNWC data files.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definitions\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 for typing information\nfrom typing import Tuple\n# Python package for time calculations\nimport pandas as pd\n# Python package for numerical calculations\nimport numpy as np\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 colorbars\nfrom matplotlib.colors import BoundaryNorm, ListedColormap\n\n# swirlspy h8 parser function\nimport swirlspy.sat.safnwc as SAFNWC\n# directory constants\nfrom swirlspy.tests.samples import DATA_DIR\nfrom swirlspy.tests.outputs import OUTPUT_DIR\n\nplt.switch_backend('agg')\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 data boundary in WGS84 (latitude)\nlatitude_from = 30.\nlatitude_to = 16.\nlongitude_from = 105.\nlongitude_to = 122.\n\nextents = (\n longitude_from, longitude_to,\n latitude_from, latitude_to\n)\n\n# Define grid size, use negative value for descending range\ngrid_size = (-.025, .025)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: Define data directory\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Supply data directory.\ndata_dir = os.path.join(DATA_DIR, \"safnwc\")\n\ninitialising_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: use ct, hrw as example. Other product will use similar flow as CT.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ct_path = os.path.join(\n data_dir,\n f'S_NWC_CT_HIMA08_HKOxRSMC-NR_20211108T000000Z.nc'\n)\n\nwith SAFNWC.read_ct(\n ct_path,\n top=latitude_from,\n bottom=latitude_to,\n y_step=grid_size[0],\n left=longitude_from,\n right=longitude_to,\n x_step=grid_size[1]\n) as ds:\n ct = ds['ct']\n ct_pal = ds['ct_pal']\n ct_attrs = ds.attrs\n\nhrw_grouping = [\n (\"all\", None)\n]\nhrw_path = os.path.join(\n data_dir,\n f'S_NWC_HRW_HIMA08_HKOxRSMC-BS_20211108T000000Z.nc'\n)\nwith SAFNWC.read_hrw(\n hrw_path,\n hrw_group=hrw_grouping,\n top=latitude_from,\n bottom=latitude_to,\n y_step=grid_size[0],\n left=longitude_from,\n right=longitude_to,\n x_step=grid_size[1]) as ds:\n hrw = ds[f\"hrw_{hrw_grouping[0][0]}\"]\n hrw_pal = ds['hrw_pal']\n hrw_attrs = ds.attrs\n\n\n# CT data\nprint(ct)\n\n# HRW data\nprint(hrw)\n\nsafnwc_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating result map\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 the crs\ncrs = ccrs.PlateCarree()\n\n# Defining coastlines\nmap_shape_file = os.path.join(DATA_DIR, \"shape/rsmc\")\ncoastline = cfeature.ShapelyFeature(\n list(shpreader.Reader(map_shape_file).geometries()),\n ccrs.PlateCarree()\n)\n\n# Defining plotting function\n\ndef _format_axis_value(v: float):\n return '%g' % (v)\n\n\ndef _setup_fig(\n axis_x: list, axis_y: list,\n time_string: str, name: str,\n norm: BoundaryNorm, cmap: ListedColormap,\n ticks: list, ticks_loc: list,\n dpi: int,\n map_shape_file: str,\n ocean_color: list,\n land_color: list,\n invert_colorbar: bool = False\n) -> Tuple[plt.Figure, plt.Axes]:\n # cofig\n width, height = 1920, 1080\n\n # color mappable\n mappable = plt.cm.ScalarMappable(\n norm=norm,\n cmap=cmap\n )\n mappable.set_array([])\n\n # base_map plotting function\n def plot_base(ax: plt.Axes):\n ax.set_extent(extents, crs=ccrs.PlateCarree())\n\n # coastline and province\n coastline_province = cfeature.ShapelyFeature(\n list(shpreader.Reader(map_shape_file).geometries()),\n ccrs.PlateCarree()\n )\n # fake the ocean color\n ax.imshow(np.tile(ocean_color, [2, 2, 1]),\n origin='upper',\n transform=ccrs.PlateCarree(),\n extent=[-180, 180, -180, 180],\n zorder=-1)\n # coastline, state, color\n ax.add_feature(coastline_province,\n facecolor=land_color, edgecolor='none', zorder=0)\n # overlay coastline, state without color\n ax.add_feature(coastline_province, facecolor='none',\n edgecolor='white', linewidth=0.5)\n\n # axis_ticks plotting function\n def plot_axis_ticks(ax: plt.Axes, fontsize: float):\n xticks = np.arange(axis_x.min(), axis_x.max(), 10)\n ax.set_xticks(xticks)\n ax.set_xticklabels(\n [f'{_format_axis_value(t)}W' if t <\n 0 else f'{_format_axis_value(t)}E' for t in xticks],\n fontsize=fontsize\n )\n\n yticks = np.arange(axis_y.min(), axis_y.max(), 10)\n ax.set_yticks(yticks)\n ax.set_yticklabels(\n [f'{_format_axis_value(t)}S' if t <\n 0 else f'{_format_axis_value(t)}N' for t in yticks],\n fontsize=fontsize\n )\n\n ##############################################################################\n # Step 3: Plotting the swirls-radar-advection, nwp-bias-corrected, blended 3 hours ahead\n #\n fig: plt.Figure = plt.figure(\n figsize=(width / dpi, height / dpi)\n )\n gs = plt.GridSpec(\n 1, 1, figure=fig,\n # wspace=0, hspace=0, top=0.95, bottom=0.05, left=0.27, right=0.845\n wspace=0, hspace=0, top=0.95, bottom=0.05, left=0.18, right=0.845\n )\n ax: plt.Axes = fig.add_subplot(\n gs[0, 0],\n projection=ccrs.PlateCarree()\n )\n fig.subplots_adjust(left=.02, right=.98, bottom=.02,\n top=.98, wspace=.0, hspace=.0)\n\n # plot base map\n plot_base(ax)\n\n # axis ticks\n fontsize = height * 0.01\n plot_axis_ticks(ax, fontsize)\n\n # title\n fontsize = height * 0.01\n ax.set_title(\n time_string,\n fontsize=fontsize,\n loc='left',\n color=\"white\"\n )\n ax.set_title(\n name,\n fontsize=fontsize,\n loc='right',\n color=\"white\"\n )\n\n fontsize = height * 0.009\n cbar_ax = fig.add_axes([0.85, 0.1, 0.01, 0.8])\n cbar_ax.tick_params(labelsize=fontsize)\n cbar = fig.colorbar(mappable, cax=cbar_ax, extend='neither')\n cbar.set_ticks(ticks_loc)\n cbar.set_ticklabels(ticks)\n\n if invert_colorbar:\n cbar.ax.invert_yaxis()\n\n fig.patch.set_facecolor('#000000')\n fig.set_edgecolor('#FFFFFF')\n for x in [ax, cbar_ax]:\n x.tick_params(color='#FFFFFF', labelcolor='#FFFFFF')\n for spine in ax.spines.values():\n spine.set_edgecolor('#FFFFFF')\n cbar.outline.set_edgecolor('#FFFFFF')\n\n return fig, ax\n\n# plot CT product\n# extract color palette information\nk = 16\ncolors = (ct_pal.data[1:k, :]) / 255.0\nlevels = np.linspace(1, k, k)\nlevel_ticks = [\n \"Cloud-free land\", \"Cloud-free sea\", \"Snow over land\", \"Sea ice\", \"Very low clouds\",\n \"Low clouds\", \"Mid-level clouds\", \"High opaque clouds\", \"Very high opaque clouds\", \"Fractional clouds\",\n \"High semitransparent thin\", \"High semi- meanly thick\", \"High semi- thick\", \"High semi- above low/med\", \"High semi- above snow/ice\"\n]\nlevel_ticks = [f'{i + 1} {l}' for i, l in enumerate(level_ticks)]\nlevel_ticks_loc = [\n l + (levels[i + 1] - l if i < k - 1 else l - levels[i - 1]) / 2\n for i, l in enumerate(levels)\n]\n\n# color map\ncmap = ListedColormap(colors)\n\n# boundary\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\n# axis\ndim_y, dim_x = ct.dims[-2:]\naxis_x = ct.coords[dim_x].values\naxis_y = ct.coords[dim_y].values\n\n# setup base figure\nfig, ax = _setup_fig(\n axis_x, axis_y,\n ct_attrs['nominal_product_time'],\n ct_attrs['long_name'],\n norm, cmap,\n level_ticks, level_ticks_loc,\n 200,\n map_shape_file,\n ocean_color=colors[1],\n land_color=colors[0]\n)\n\n# plot actual data on map\nct.where(ct >= 2).plot(\n ax=ax,\n cmap=cmap,\n norm=norm,\n add_colorbar=False,\n add_labels=False\n)\n\nfig.savefig(\n os.path.join(OUTPUT_DIR, \"safnwc_ct.png\"),\n dpi=200,\n bbox_inches=\"tight\",\n pad_inches=0.1,\n facecolor=fig.get_facecolor(),\n edgecolor='none'\n)\n\n# plot HRW product\n\n# axis\ndim_y, dim_x = hrw.dims[-2:]\naxis_x = hrw.coords[dim_x].values\naxis_y = hrw.coords[dim_y].values\n\n# extract necessary data\nbarb = list(hrw.coords['barb'].values)\nap_idx = barb.index('air_pressure')\nu_idx = barb.index('u')\nv_idx = barb.index('v')\n\n# filter non null data\ny = np.repeat(axis_y, hrw.shape[2]).reshape(hrw.shape[1:])\nx = np.tile(axis_x, hrw.shape[1]).reshape(hrw.shape[1:])\nidx = np.isnan(hrw[ap_idx, :, :])\nhrw_data = np.concatenate((hrw, [y], [x]))[:, ~idx]\n\n# filter non null data\nx = hrw_data[-1, :]\ny = hrw_data[-2, :]\nu = hrw_data[u_idx, :]\nv = hrw_data[v_idx, :]\nflip_barb = y < 0\nair_pressure = hrw_data[ap_idx, :] / 100\n\n# extract color palette information\nlevels = [150, 200, 250, 300, 350, 400, 450, 500,\n 550, 600, 650, 700, 750, 800, 850, 900, 950]\ncolors = (hrw_pal.data[11:199:11]) / 255.0\n\n# color map\ncmap = ListedColormap(colors)\n# boundary\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\nfig, ax = _setup_fig(\n axis_x, axis_y,\n hrw_attrs['nominal_product_time'],\n f\"High Resolution Winds: 100-1000 hPa\",\n norm, cmap,\n levels, levels,\n 200,\n map_shape_file,\n ocean_color=np.array([178, 208, 254], dtype=np.uint8) / 255.,\n land_color=cfeature.COLORS['land'],\n invert_colorbar=True\n)\n\nax = ax.barbs(\n x, y, u, v, air_pressure,\n cmap=cmap, norm=norm,\n sizes={'spacing': 0.2, 'width': 0.15, 'height': 0.3, 'emptybarb': 0.0},\n length=3.5, linewidth=0.25,\n flip_barb=flip_barb,\n pivot='middle',\n barb_increments={'half': 2.5722810989, 'full': 5.1445621978, 'flag': 25.7228109888}\n)\n\nfig.savefig(\n os.path.join(OUTPUT_DIR, \"safnwc_hrw.png\"),\n dpi=200,\n bbox_inches=\"tight\",\n pad_inches=0.1,\n facecolor=fig.get_facecolor(),\n edgecolor='none'\n)\n\nsafnwc_image_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\"SAFNWC data parsing time: {safnwc_time}\")\nprint(f\"Plotting safnwc image time: {safnwc_image_time}\")\n\nprint(f\"Time to initialise: {initialising_time - start_time}\")\nprint(f\"Time to run data parsing: {safnwc_time - initialising_time}\")\nprint(f\"Time to plot data images: {safnwc_image_time - safnwc_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 }