examples/python/Accumulated_prec_amount_from_radar_GeoTiff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate precipitation sum from GeoTiff radar images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finnish Meteorological Institute is providing radar images as GeoTiff in S3 bucket. The data is free and openly available with CC4BY license.\n",
    "\n",
    "This is a short and simple example how to calculate accumulated precipitation amount from GeoTiff radar images of precipitation rate (rr) with 5 minute interval. The example is ment to illustrate how to use the data with python.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we install necessary libraries. Rasterio (https://github.com/mapbox/rasterio) is used to fetch the data and pyproj (https://pypi.org/project/pyproj/) to calculate coordinate transformations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install rasterio[s3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio as rio\n",
    "from datetime import datetime, timedelta\n",
    "from pyproj import Proj, Transformer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Small helper function to handle time steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roundTime(dt, roundTo=300):\n",
    "   \"\"\"\n",
    "   Round a datetime object to any time lapse in seconds\n",
    "   \n",
    "   dt      : datetime.datetime object\n",
    "             DateTime to round\n",
    "   roundTo : int\n",
    "             Closest number of seconds to round to, default 1 minute. \n",
    "   \n",
    "   Returns\n",
    "   -------\n",
    "   datetime.datetime object\n",
    "            Modified datetime   \n",
    "   \"\"\"\n",
    "   seconds = (dt.replace(tzinfo=None) - dt.min).seconds\n",
    "   rounding = (seconds+roundTo/2) // roundTo * roundTo\n",
    "   return dt + timedelta(0,rounding-seconds,-dt.microsecond)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following function opens the geotiff image from S3 bucket and sample pixel value from requested coordinates. Note that images are in EPSG:3067 projection so we need to reproject the coordinates. \n",
    "\n",
    "Finally, pixel values are converted to millimeters. Conversions are documented here: https://en.ilmatieteenlaitos.fi/open-data-manual-radar-data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_timestep(bucket, key, x, y, mode='rr', timestep=300):\n",
    "    \"\"\" Process time step\n",
    "    \n",
    "    bucket   : str \n",
    "               S3 bucket name (where radar images exist)\n",
    "    key      : str\n",
    "               Radar image key in the bucket\n",
    "    lat      : float\n",
    "               Latitude of the point of interest\n",
    "    lon      : float\n",
    "               Longitude of the point of interest    \n",
    "    mode     : str\n",
    "               Defines how pixel values are converted to millimeters\n",
    "    timestep : int\n",
    "               Time step of radar images (used in pixel --> mm conversion)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float \n",
    "         Precipitation amount in millimeters\n",
    "    \"\"\"    \n",
    "    with rio.open('s3://{}/{}'.format(bucket, key)) as src:\n",
    "        for val in src.sample([(x, y)]):\n",
    "            if mode == 'rr':\n",
    "                # From https://en.ilmatieteenlaitos.fi/open-data-manual-radar-data\n",
    "                return val[0]*.01*timestep/3600\n",
    "            elif mode == 'dbz':\n",
    "                # From https://en.wikipedia.org/wiki/DBZ_(meteorology)\n",
    "                dbz = 0.5 * val[0] - 32\n",
    "                return (10**(dbz/10)/200)**(5/8)*timestep/3600\n",
    "            else:\n",
    "                raise Exception('mode not implemented')\n",
    "        \n",
    "    return 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last, we define a function to go through requested time range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_time_range(bucket, starttime, endtime, lat, lon, timestep, mode = 'rr', filename_start='FIN-DBZ-3067-250M/', filename_end='_FIN-DBZ-3067-250M.tif'):\n",
    "    \"\"\" Process time range \n",
    "    bucket    : str \n",
    "                S3 bucket name (where radar images exist)\n",
    "    starttime : str\n",
    "                Start time of the time range in format %Y-%m-%d %H:%M:%S\n",
    "    endtime   : str\n",
    "                End time of the time range in format %Y-%m-%d %H:%M:%S                \n",
    "    lat       : float\n",
    "                Latitude of the point of interest\n",
    "    lon       : float\n",
    "                Longitude of the point of interest    \n",
    "    timestep  : int\n",
    "                Time step of radar images\n",
    "    mode      : str\n",
    "                Defines how pixel values are converted to millimeters\n",
    "    filename_start : str\n",
    "                     For example path of the image. Images are assumed to found with name \n",
    "                     filename_start+timestamp+filename_end.  \n",
    "    filename_start : str\n",
    "                 Images are assumed to found with name filename_start+timestamp+filename_end.  \n",
    "                      \n",
    "    Returns\n",
    "    -------\n",
    "    float \n",
    "         Precipitation amount in millimeters\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    start = roundTime(datetime.strptime(starttime, '%Y-%m-%d %H:%M:%S'))\n",
    "    end = roundTime(datetime.strptime(endtime, '%Y-%m-%d %H:%M:%S'))    \n",
    "    time_it = start    \n",
    "    \n",
    "    transformer = Transformer.from_crs(\"epsg:4326\", \"epsg:3067\")\n",
    "    x,y = transformer.transform(lat, lon)\n",
    "    \n",
    "    prec_sum = 0\n",
    "    while time_it <= end:\n",
    "        key = '{year}/{month}/{day}/{prefix}{stamp}{suffix}'.format(year=time_it.strftime('%Y'),\n",
    "                                                                    month=time_it.strftime('%m'),\n",
    "                                                                    day=time_it.strftime('%d'),\n",
    "                                                                    stamp=time_it.strftime('%Y%m%d%H%M'),\n",
    "                                                                    prefix=filename_start,\n",
    "                                                                    suffix=filename_end)\n",
    "        prec_sum += process_timestep(bucket, key, x, y, timestep=timestep)\n",
    "\n",
    "        if time_it.minute == 0: print('{:.2f}'.format(prec_sum), end='')\n",
    "        else: print('.', end='')\n",
    "            \n",
    "        time_it +=  timedelta(seconds=timestep)\n",
    "\n",
    "    print('')\n",
    "    return prec_sum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we run everything:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "bucket = 'fmi-opendata-radar-geotiff'\n",
    "starttime = '2020-09-06 00:00:00'\n",
    "endtime = '2020-09-06 23:59:00'\n",
    "lat = 63.1593\n",
    "lon = 29.8346"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.06...........0.10...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.17...........0.21...........0.24...........0.56...........0.91...........1.32...........1.67\n",
      "Precipitation sum at 63.1593,29.8346 (Koli, Finland) on 11th July 2020 is 1.67 mm\n",
      "CPU times: user 8.6 s, sys: 389 ms, total: 8.99 s\n",
      "Wall time: 1min 53s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "prec_sum = process_time_range(bucket, starttime, endtime, lat, lon, 300, 'dbz')\n",
    "print('Precipitation sum at {},{} (Koli, Finland) on 11th July 2020 is {:.2f} mm'.format(lat,lon,prec_sum))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dask version\n",
    "\n",
    "As the above computation contains almost 300 geotiff images to process, it does take some time. If we want to calculate precipitation sum over longer period or do some other heavy computation, it may be good idea to distritbute the calculation. Following example shows how this can achieved using Dask (https://docs.dask.org/en/latest/). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to install the envinronment. This example does not cover installing dask but all assets used to setup the dask cluster are included in the repository. A good tutorial for setting up the cluster can be found at: https://towardsdatascience.com/serverless-distributed-data-pre-processing-using-dask-amazon-ecs-and-python-part-1-a6108c728cc4\n",
    "\n",
    "To use the cluster, we need several libraries to be installed. Following commands install the require libraries. Versions of each libraries are defined as they need to mach with the cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip uninstall -y enum34\n",
    "!pip install blosc==1.8.3\n",
    "!pip install --upgrade joblib==0.16 dask[dataframe]==2.25.0 distributed==2.25.0 lz4==3.1.0 cloudpickle==1.6.0 msgpack numpy==1.19.1 toolz tornado blosc==1.7.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we import necessary libraries and initialize dask client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client\n",
    "from dask import delayed\n",
    "import dask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border: 2px solid white;\">\n",
       "<tr>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3 style=\"text-align: left;\">Client</h3>\n",
       "<ul style=\"text-align: left; list-style: none; margin: 0; padding: 0;\">\n",
       "  <li><b>Scheduler: </b>tcp://Dask-Scheduler.local-dask:8786</li>\n",
       "  <li><b>Dashboard: </b><a href='http://Dask-Scheduler.local-dask:8787/status' target='_blank'>http://Dask-Scheduler.local-dask:8787/status</a></li>\n",
       "</ul>\n",
       "</td>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3 style=\"text-align: left;\">Cluster</h3>\n",
       "<ul style=\"text-align: left; list-style:none; margin: 0; padding: 0;\">\n",
       "  <li><b>Workers: </b>2</li>\n",
       "  <li><b>Cores: </b>8</li>\n",
       "  <li><b>Memory: </b>34.36 GB</li>\n",
       "</ul>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Client: 'tcp://172.31.35.74:8786' processes=2 threads=8, memory=34.36 GB>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "client = Client('Dask-Scheduler.local-dask:8786')\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define necessary variables. Here we use the same values with previous example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "bucket = 'fmi-opendata-radar-geotiff'\n",
    "lat = 63.1593\n",
    "lon = 29.8346\n",
    "timestep=300\n",
    "filename_start='FIN-DBZ-3067-250M/'\n",
    "filename_end='_FIN-DBZ-3067-250M.tif'\n",
    "starttime = datetime.strptime('2020-09-06 00:00:00', \"%Y-%m-%d %H:%M:%S\")\n",
    "endtime = datetime.strptime('2020-09-06 23:59:00', \"%Y-%m-%d %H:%M:%S\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing is also almost similar to local version. We go through the time steps but instead of directly calculating the value, we define the calculation and call dask to compute all at once. Dask handles the distribution and gathering results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Precipitation sum at 63.1593,29.8346 (Koli, Finland) on between 06th Sep 2020 - 06th Sep 2020 is 0.08 mm\n",
      "CPU times: user 693 ms, sys: 33.1 ms, total: 726 ms\n",
      "Wall time: 28.8 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "transformer = Transformer.from_crs(\"epsg:4326\", \"epsg:3067\")\n",
    "x,y = transformer.transform(lat, lon)\n",
    "\n",
    "dfs, df = [], []\n",
    "time_it = starttime\n",
    "while time_it <= endtime:\n",
    "    key = '{year}/{month}/{day}/{prefix}{stamp}{suffix}'.format(year=time_it.strftime('%Y'),\n",
    "                                                                month=time_it.strftime('%m'),\n",
    "                                                                day=time_it.strftime('%d'),\n",
    "                                                                stamp=time_it.strftime('%Y%m%d%H%M'),\n",
    "                                                                prefix=filename_start,\n",
    "                                                                suffix=filename_end)\n",
    "    dfs.append(delayed(process_timestep)(bucket, key, x, y, timestep=timestep, mode='dbz'))\n",
    "    time_it +=  timedelta(seconds=timestep)\n",
    "\n",
    "df = dask.compute(*dfs)\n",
    "\n",
    "print('Precipitation sum at {},{} (Koli, Finland) on between {} - {} is {:.2f} mm'.format(lat,\n",
    "                                                                                          lon,\n",
    "                                                                                          starttime.strftime('%dth %b %Y'), \n",
    "                                                                                          endtime.strftime('%dth %b %Y'), \n",
    "                                                                                          sum(df)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "conda_python3",
   "language": "python",
   "name": "conda_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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}