{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Wind energy production forecast from Met.no weather forecast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this notebook we illustrate how one might improve weather forecast for wind energy production, considering that the height of wind turbines doesn’t match the height of wind speed commonly used in weather forecasts. Also, that wind energy production does not depend only on wind speed, but also air density.\n",
    "\n",
    "In this Notebook we will use three weather forecast from __[Planet OS Datahub](data.planetos.com)__:\n",
    "\n",
    "1. __[GFS](http://data.planetos.com/datasets/noaa_gfs_pgrb2_global_forecast_recompute_0.25degree)__ - A global model produced by NOAA that provides a 10-day forecast at a resolution of 0.25 degrees (15 days with lower resolution and fewer variables);\n",
    "\n",
    "2. __[FMI HIRLAM](http://data.planetos.com/datasets/fmi_hirlam_surface)__ - A regional model limited to Northern-Europe and Northern-Atlantic with a forecast period of two days and resolution of 0.07 degrees;\n",
    "\n",
    "3. __[Met.no HARMONIE](http://data.planetos.com/datasets/metno_harmonie_metcoop)__ & __[Met.no HARMONIE dedicated wind forecast](http://data.planetos.com/datasets/metno_harmonie_wind_det)__ A regional model limited to Scandinavia with a resolution of 0.025 x 0.05 degrees. The dedicated wind forecast version is augmented by Planet OS to include the air density as explicit variable.\n",
    "\n",
    "\n",
    "_API documentation is available at http://docs.planetos.com. If you have questions or comments, join the Planet OS Slack community to chat with our development team. For general information on usage of IPython/Jupyter and Matplotlib, please refer to their corresponding documentation. https://ipython.org/ and http://matplotlib.org/_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import urllib.request\n",
    "import numpy as np\n",
    "import simplejson as json\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "import warnings\n",
    "import datetime\n",
    "import dateutil.parser\n",
    "import matplotlib.cbook\n",
    "warnings.filterwarnings(\"ignore\",category=matplotlib.cbook.mplDeprecation)\n",
    "import requests\n",
    "from netCDF4 import Dataset\n",
    "from dh_py_access.lib.dataset import dataset as dataset\n",
    "import dh_py_access.package_api as package_api\n",
    "import dh_py_access.lib.datahub as datahub\n",
    "np.warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>Please put your datahub API key into a file called APIKEY and place it to the notebook folder or assign your API key directly to the variable API_key!</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "API_key = open(\"APIKEY\").read().strip()\n",
    "dh=datahub.datahub_main(API_key)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a simple dh_py_access library for controlling REST requests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fmi_hirlam_surface = dataset('fmi_hirlam_surface',dh)\n",
    "metno_harmonie_metcoop = dataset('metno_harmonie_metcoop',dh)\n",
    "metno_harmonie_wind = dataset('metno_harmonie_wind_det',dh)\n",
    "gfs = dataset('noaa_gfs_pgrb2_global_forecast_recompute_0.25degree',dh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "dh_py_access provides some easy to use functions to list dataset metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Does not look good in github \n",
    "##gfs.variable_names()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize coordinates and reftime. Please note that more recent reftime, than computed below, may be available for particular dataset. We just use conservative example for demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon = 22\n",
    "lat = 59+6./60\n",
    "today = datetime.datetime.today()\n",
    "reft = datetime.datetime(today.year,today.month,today.day,int(today.hour/6)*6) - datetime.timedelta(hours=12)\n",
    "reft = reft.isoformat()\n",
    "##reft = \"2018-02-11T18:00:00\"\n",
    "arg_dict = {'lon':lon,'lat':lat,'reftime_start':reft,'reftime_end':reft,'count':250}\n",
    "arg_dict_metno_wind_det = dict(arg_dict, **{'vars':'wind_u_z,wind_v_z,air_density_z'})\n",
    "arg_dict_metno_harm_metcoop = dict(arg_dict, **{'vars':'u_wind_10m,v_wind_10m'})\n",
    "arg_dict_hirlam = dict(arg_dict, **{'vars':'u-component_of_wind_height_above_ground,v-component_of_wind_height_above_ground'})\n",
    "arg_dict_gfs = dict(arg_dict, **{'vars':'ugrd_m,vgrd_m','count':450})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fetch data and convert to Pandas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "\"['context' 'axes'] not in index\"",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-345fcf5c9588>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdmw\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmetno_harmonie_wind\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_json_data_in_pandas\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0marg_dict_metno_wind_det\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mdmm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmetno_harmonie_metcoop\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_json_data_in_pandas\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0marg_dict_metno_harm_metcoop\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mdhs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfmi_hirlam_surface\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_json_data_in_pandas\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0marg_dict_hirlam\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mdgfs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgfs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_json_data_in_pandas\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0marg_dict_gfs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/etoodu/Desktop/planetOS/git/notebooks/api-examples/dh_py_access/lib/dataset.py\u001b[0m in \u001b[0;36mget_json_data_in_pandas\u001b[0;34m(self, count, z, pandas, **kwargs)\u001b[0m\n\u001b[1;32m    103\u001b[0m             \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"z\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    104\u001b[0m         \u001b[0mretjson\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mparse_urls\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdatahub\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mserver\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdatahub\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mversion\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"datasets/{0}/point\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdatasetkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdatahub\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapikey\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mclean_reftime\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjson\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 105\u001b[0;31m         \u001b[0;32mif\u001b[0m \u001b[0mpandas\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mretjson\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mconvert_json_to_some_pandas\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mretjson\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'entries'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    106\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mretjson\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    107\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/etoodu/Desktop/planetOS/git/notebooks/api-examples/dh_py_access/lib/dataset.py\u001b[0m in \u001b[0;36mconvert_json_to_some_pandas\u001b[0;34m(injson)\u001b[0m\n\u001b[1;32m     95\u001b[0m             \u001b[0;34m[\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnew_dict\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'axes'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'axes'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnew_dict\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'data'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'data'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0minjson\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     96\u001b[0m             \u001b[0mpd_temp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minjson\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 97\u001b[0;31m             \u001b[0mdev_frame\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd_temp\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'context'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'axes'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnew_dict\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mparam_list\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     98\u001b[0m             \u001b[0mdev_frame\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdev_frame\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdev_frame\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'reftime'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mdev_frame\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'reftime'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     99\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mdev_frame\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/etoodu/anaconda/envs/py36/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   1956\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mSeries\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mIndex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1957\u001b[0m             \u001b[0;31m# either boolean or fancy integer index\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1958\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_array\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1959\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mDataFrame\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1960\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_frame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/etoodu/anaconda/envs/py36/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m_getitem_array\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   2000\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtake\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconvert\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2001\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2002\u001b[0;31m             \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_convert_to_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2003\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtake\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconvert\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2004\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/etoodu/anaconda/envs/py36/lib/python3.6/site-packages/pandas/core/indexing.py\u001b[0m in \u001b[0;36m_convert_to_indexer\u001b[0;34m(self, obj, axis, is_setter)\u001b[0m\n\u001b[1;32m   1229\u001b[0m                 \u001b[0mmask\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcheck\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1230\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mmask\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0many\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1231\u001b[0;31m                     \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'%s not in index'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mobjarr\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mmask\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1232\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1233\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0m_values_from_object\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: \"['context' 'axes'] not in index\""
     ]
    }
   ],
   "source": [
    "dmw = metno_harmonie_wind.get_json_data_in_pandas(**arg_dict_metno_wind_det)\n",
    "dmm = metno_harmonie_metcoop.get_json_data_in_pandas(**arg_dict_metno_harm_metcoop)\n",
    "dhs = fmi_hirlam_surface.get_json_data_in_pandas(**arg_dict_hirlam)\n",
    "dgfs = gfs.get_json_data_in_pandas(**arg_dict_gfs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## show how to filter Pandas\n",
    "## dgfs[dgfs['z']==80]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter out necessary data from the DataFrames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vel80_metno = np.array(np.sqrt(dmw[dmw['z']==80]['wind_u_z']**2 + dmw[dmw['z']==80]['wind_v_z']**2))\n",
    "vel10_metno = np.array(np.sqrt(dmm['u_wind_10m']**2 + dmm['v_wind_10m']**2))\n",
    "vel10_hirlam = np.array(np.sqrt(dhs['u-component_of_wind_height_above_ground']**2 + \n",
    "                          dhs['v-component_of_wind_height_above_ground']**2))\n",
    "vel10_gfs = np.sqrt(dgfs[dgfs['z']==10]['ugrd_m']**2+dgfs[dgfs['z']==10]['vgrd_m']**2)\n",
    "vel80_gfs = np.sqrt(dgfs[dgfs['z']==80]['ugrd_m']**2+dgfs[dgfs['z']==80]['vgrd_m']**2)\n",
    "t_metno = [dateutil.parser.parse(i) for i in dmw[dmw['z']==80]['time']]\n",
    "t_metno_10 = [dateutil.parser.parse(i) for i in dmm['time']]\n",
    "t_hirlam = [dateutil.parser.parse(i) for i in dhs['time']]\n",
    "t_gfs_10 = [dateutil.parser.parse(i) for i in dgfs[dgfs['z']==10]['time']]\n",
    "t_gfs_80 = [dateutil.parser.parse(i) for i in dgfs[dgfs['z']==80]['time']]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "days = mdates.DayLocator()\n",
    "daysFmt = mdates.DateFormatter('%Y-%m-%d')\n",
    "hours = mdates.HourLocator()\n",
    "ax.set_ylabel(\"wind speed\")\n",
    "ax.plot(t_metno, vel80_metno, label='Metno 80m')\n",
    "ax.plot(t_metno_10, vel10_metno, label='Metno 10m')\n",
    "ax.plot(t_hirlam, vel10_hirlam, label='HIRLAM 10m')\n",
    "gfs_lim=67\n",
    "ax.plot(t_gfs_10[:gfs_lim], vel10_gfs[:gfs_lim], label='GFS 10m')\n",
    "ax.plot(t_gfs_80[:gfs_lim], vel80_gfs[:gfs_lim], label='GFS 80m')\n",
    "ax.xaxis.set_major_locator(days)\n",
    "ax.xaxis.set_major_formatter(daysFmt)\n",
    "ax.xaxis.set_minor_locator(hours)\n",
    "#fig.autofmt_xdate()\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.savefig(\"model_comp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We easily see that differences between models can be larger than difference of 10m to 80m winds in the same model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What role does density play in energy production?\n",
    "From the theory we know that wind energy production is roughly \n",
    "\n",
    "$1/2 A \\rho \\mathbf{v}^3$, \n",
    "\n",
    "where $A$ is area, $\\rho$ is air density and $\\mathbf{v}$ is wind speed. We are not concerned about $A$, which is a turbine parameter, but we can analyse the linear relation of density and cube relation of wind speed itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's see how the density varies over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax2 = ax.twinx()\n",
    "\n",
    "\n",
    "days = mdates.DayLocator()\n",
    "daysFmt = mdates.DateFormatter('%Y-%m-%d')\n",
    "hours = mdates.HourLocator()\n",
    "ax.plot(t_metno,vel80_metno)\n",
    "\n",
    "aird80 = dmw[dmw['z']==80]['air_density_z']\n",
    "ax2.plot(t_metno,aird80,c='g')\n",
    "ax.xaxis.set_major_locator(days)\n",
    "ax.xaxis.set_major_formatter(daysFmt)\n",
    "ax.xaxis.set_minor_locator(hours)\n",
    "ax.set_ylabel(\"wind speed\")\n",
    "ax2.set_ylabel(\"air density\")\n",
    "fig.tight_layout()\n",
    "#fig.autofmt_xdate()\n",
    "plt.savefig(\"density_80m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No let's energy production looks compared to the wind speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax2 = ax.twinx()\n",
    "ax2.set_ylabel(\"energy production\")\n",
    "ax.set_ylabel(\"wind speed\")\n",
    "ax.plot(t_metno,vel80_metno, c='b', label='wind speed')\n",
    "ax2.plot(t_metno,aird80*vel80_metno**3, c='r', label='energy prod w.dens')\n",
    "ax.xaxis.set_major_locator(days)\n",
    "ax.xaxis.set_major_formatter(daysFmt)\n",
    "ax.xaxis.set_minor_locator(hours)\n",
    "fig.autofmt_xdate()\n",
    "lines, labels = ax.get_legend_handles_labels()\n",
    "lines2, labels2 = ax2.get_legend_handles_labels()\n",
    "ax2.legend(lines + lines2, labels + labels2, loc=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's analyse how much density changes vary over the whole domain during one forecast. For this purpose, we download the whole density field with package api"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "density = package_api.package_api(dh,'metno_harmonie_wind_det','air_density_z',-20,60,10,80,'full_domain_harmonie') \n",
    "density.make_package()\n",
    "density.download_package()\n",
    "density_data = Dataset(density.get_local_file_name())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## biggest change of density in one location during forecast period\n",
    "maxval = np.nanmax(density_data.variables['air_density_z'],axis=0)\n",
    "minval = np.nanmin(density_data.variables['air_density_z'],axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maximum relative change of air density in single location is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(np.nanmax(maxval-minval),np.nanmean(maxval-minval))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}