{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Linear blending\n\nThis tutorial shows how to construct a simple linear blending between a STEPS\nensemble nowcast and a Numerical Weather Prediction (NWP) rainfall forecast. The\nused datasets are from the Bureau of Meteorology, Australia.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nfrom datetime import datetime\n\nfrom matplotlib import pyplot as plt\n\nimport pysteps\nfrom pysteps import io, rcparams, blending\nfrom pysteps.visualization import plot_precip_field" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read the radar images and the NWP forecast\n\nFirst, we import a sequence of 3 images of 10-minute radar composites\nand the corresponding NWP rainfall forecast that was available at that time.\n\nYou need the pysteps-data archive downloaded and the pystepsrc file\nconfigured with the data_source paths pointing to data folders.\nAdditionally, the pysteps-nwp-importers plugin needs to be installed, see\nhttps://github.com/pySTEPS/pysteps-nwp-importers.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Selected case\ndate_radar = datetime.strptime(\"202010310400\", \"%Y%m%d%H%M\")\n# The last NWP forecast was issued at 00:00\ndate_nwp = datetime.strptime(\"202010310000\", \"%Y%m%d%H%M\")\nradar_data_source = rcparams.data_sources[\"bom\"]\nnwp_data_source = rcparams.data_sources[\"bom_nwp\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load the data from the archive\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "root_path = radar_data_source[\"root_path\"]\npath_fmt = \"prcp-c10/66/%Y/%m/%d\"\nfn_pattern = \"66_%Y%m%d_%H%M00.prcp-c10\"\nfn_ext = radar_data_source[\"fn_ext\"]\nimporter_name = radar_data_source[\"importer\"]\nimporter_kwargs = radar_data_source[\"importer_kwargs\"]\ntimestep = 10.0\n\n# Find the radar files in the archive\nfns = io.find_by_date(\n date_radar, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2\n)\n\n# Read the radar composites\nimporter = io.get_method(importer_name, \"importer\")\nradar_precip, _, radar_metadata = io.read_timeseries(fns, importer, **importer_kwargs)\n\n# Import the NWP data\nfilename = os.path.join(\n nwp_data_source[\"root_path\"],\n datetime.strftime(date_nwp, nwp_data_source[\"path_fmt\"]),\n datetime.strftime(date_nwp, nwp_data_source[\"fn_pattern\"])\n + \".\"\n + nwp_data_source[\"fn_ext\"],\n)\n\nnwp_importer = io.get_method(\"bom_nwp\", \"importer\")\nnwp_precip, _, nwp_metadata = nwp_importer(filename)\n\n# Only keep the NWP forecasts from the last radar observation time (2020-10-31 04:00)\n# End of the forecast is 18 time steps (+3 hours) in advance.\nprecip_nwp = nwp_precip[24:43, :, :]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Pre-processing steps\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Make sure the units are in mm/h\nconverter = pysteps.utils.get_method(\"mm/h\")\nradar_precip, radar_metadata = converter(radar_precip, radar_metadata)\nprecip_nwp, nwp_metadata = converter(precip_nwp, nwp_metadata)\n\n# Threshold the data\nradar_precip[radar_precip < 0.1] = 0.0\nprecip_nwp[precip_nwp < 0.1] = 0.0\n\n# Plot the radar rainfall field and the first time step of the NWP forecast.\n# For the initial time step (t=0), the NWP rainfall forecast is not that different\n# from the observed radar rainfall, but it misses some of the locations and\n# shapes of the observed rainfall fields. Therefore, the NWP rainfall forecast will\n# initially get a low weight in the blending process.\ndate_str = datetime.strftime(date_radar, \"%Y-%m-%d %H:%M\")\nplt.figure(figsize=(10, 5))\nplt.subplot(121)\nplot_precip_field(\n radar_precip[-1, :, :],\n geodata=radar_metadata,\n title=f\"Radar observation at {date_str}\",\n)\nplt.subplot(122)\nplot_precip_field(\n precip_nwp[0, :, :], geodata=nwp_metadata, title=f\"NWP forecast at {date_str}\"\n)\nplt.tight_layout()\nplt.show()\n\n# Only keep the NWP forecasts from 2020-10-31 04:05 onwards, because the first\n# forecast lead time starts at 04:05.\nprecip_nwp = precip_nwp[1:]\n\n# Transform the radar data to dB - this transformation is useful for the motion\n# field estimation and the subsequent nowcasts. The NWP forecast is not\n# transformed, because the linear blending code sets everything back in mm/h\n# after the nowcast.\ntransformer = pysteps.utils.get_method(\"dB\")\nradar_precip, radar_metadata = transformer(radar_precip, radar_metadata, threshold=0.1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Determine the velocity field for the radar rainfall nowcast\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "oflow_method = pysteps.motion.get_method(\"lucaskanade\")\nvelocity_radar = oflow_method(radar_precip)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The linear blending of nowcast and NWP rainfall forecast\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Calculate the blended precipitation field\nprecip_blended = blending.linear_blending.forecast(\n precip=radar_precip[-1, :, :],\n precip_metadata=radar_metadata,\n velocity=velocity_radar,\n timesteps=18,\n timestep=10,\n nowcast_method=\"extrapolation\", # simple advection nowcast\n precip_nwp=precip_nwp,\n precip_nwp_metadata=nwp_metadata,\n start_blending=60, # in minutes (this is an arbritrary choice)\n end_blending=120, # in minutes (this is an arbritrary choice)\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualize the output\n\nThe linear blending starts at 60 min, so during the first 60 minutes the\nblended forecast only consists of the extrapolation forecast (consisting of an\nextrapolation nowcast). Between 60 and 120 min, the NWP forecast gradually gets more\nweight, whereas the extrapolation forecasts gradually gets less weight.\nAfter 120 min, the blended forecast entirely consists of the NWP rainfall\nforecast.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = plt.figure(figsize=(4, 8))\n\nleadtimes_min = [30, 60, 90, 120]\nn_leadtimes = len(leadtimes_min)\nfor n, leadtime in enumerate(leadtimes_min):\n\n # Nowcast with blending into NWP\n plt.subplot(n_leadtimes, 2, n * 2 + 1)\n plot_precip_field(\n precip_blended[int(leadtime / timestep) - 1, :, :],\n geodata=radar_metadata,\n title=f\"Nowcast +{leadtime} min\",\n axis=\"off\",\n colorbar=False,\n )\n\n # Raw NWP forecast\n plt.subplot(n_leadtimes, 2, n * 2 + 2)\n plot_precip_field(\n precip_nwp[int(leadtime / timestep) - 1, :, :],\n geodata=nwp_metadata,\n title=f\"NWP +{leadtime} min\",\n axis=\"off\",\n colorbar=False,\n )\n\nplt.tight_layout()\nplt.show()" ] } ], "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.8.6" } }, "nbformat": 4, "nbformat_minor": 0 }