{
 "metadata": {
  "name": "delwaq"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "\n",
      "# ---------------------------------------------------------------------------\n",
      "# Created by Gerrit Hendriksen (gerrit.hendriksen@deltares.nl)\n",
      "# v1.0 created on 10-03-2012 (ddmmyyyy)\n",
      "# v1.1 adapted on 15-03-2012 (ddmmyyyy) added matplotlib to module import for\n",
      "#      display of graphs\n",
      "#\n",
      "# Description: Open netCDF4 on OPeNDAP for reading, querying and display time\n",
      "#              series on via matplotlib\n",
      "# ---------------------------------------------------------------------------\n",
      "\n",
      "# load necessary libs\n",
      "\n",
      "from multiprocessing import Pool\n",
      "\n",
      "import numpy as np\n",
      "\n",
      "import netCDF4\n",
      "import netcdftime\n",
      "\n",
      "import pandas\n",
      "\n",
      "from shapely.geometry import Polygon,Point\n",
      "\n",
      "import openearthtools.io.opendap.catalog"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Users/fedorbaart/.virtualenvs/main/lib/python2.7/site-packages/pytz/__init__.py:35: UserWarning: Module argparse was already imported from /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/argparse.pyc, but /Users/fedorbaart/.virtualenvs/main/lib/python2.7/site-packages is being added to sys.path\n",
        "  from pkg_resources import resource_stream\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def mean_concentration(url_indices):\n",
      "    \"\"\"compute the mean concentration for a url\"\"\"\n",
      "    url, indices = url_indices\n",
      "    dataset = netCDF4.Dataset(url)\n",
      "    times =  netcdftime.num2date(dataset.variables['time'], dataset.variables['time'].units)\n",
      "    conc = dataset.variables[\"concentration\"][:]\n",
      "    rows = []\n",
      "    for tidx, t in enumerate(times):\n",
      "        if indices.shape[0] < 1:\n",
      "            continue\n",
      "        arrconc = conc[tidx,indices[:,0], indices[:,1]]\n",
      "        rows.append((t, t.strftime('%Y-%m'), arrconc.mean()))\n",
      "    dataset.close()\n",
      "    return rows"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def delwaq_species():\n",
      "    # Namen van Delwaq parameters ter vergelijking:\n",
      "    #Delwaq name:Long name,Unit,CF-name\n",
      "    dctspecies = {'NH4':['Ammonium (NH4)','gN/m3','mass_concentration_of_ammonium_in_sea_water'],\n",
      "                  'NO3':['Nitrate (NO3)','gN/m3','mass_concentration_of_nitrate_in_sea_water'],\n",
      "                  'PO4':['Ortho-Phosphate (PO4)','gP/m3','mass_concentration_of_phosphate_in_sea_water'],\n",
      "                  'Si':['dissolved Silica (Si)','gSi/m3','mass_concentration_of_silicate_in_sea_water'],\n",
      "                  'ExtVl':['total extinction coefficient visible light','1/m','volume_attenuation_coefficient_of_downwelling_radiative_flux_in_sea_water'],\n",
      "                  'Temp':['ambient water temperature','degrees C','sea_water_temperature_expressed_as_degrees_celsius'],\n",
      "                  'fPPtot':['total net primary production','gC/m2/d','tendency_of_mass_concentration_of_particulate_organic_matter_expressed_as_carbon_in_sea_water_due_to_net_primary_production'],\n",
      "                  'Phyt':['total carbon in phytoplankton','gC/m3','mass_concentration_of_phytoplankton_expressed_as_carbon_in_sea_water'],\n",
      "                  'Chlfa':['Chlorophyll-a concentration','g/m3','mass_concentration_of_chlorophyll_in_sea_water'],\n",
      "                  'Chlfa2':['Chlorophyll-a concentration','g/m2','mass_of_chlorophyll_in_sea_water_per_unit_area']\n",
      "                 }\n",
      "    return dctspecies\n",
      "\n",
      "# credentials for VECTORS database (i.e. readonly)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "import urllib\n",
      "\n",
      "apoly=Polygon()\n",
      "aspecies='P04'\n",
      "\n",
      "\n",
      "\n",
      "# TODO check why this is not used\n",
      "statsq='34F3'\n",
      "\n",
      "catalogurl = 'http://opendap.deltares.nl/thredds/catalog/opendap/deltares/delwaq/catalog.xml'\n",
      "urls = list(openearthtools.io.opendap.catalog.getchildren(catalogurl))\n",
      "\n",
      "print urls[0]\n",
      "# lookup space for the first dataset\n",
      "dataset = netCDF4.Dataset(urls[0])\n",
      "lon = dataset.variables[\"lon\"][:]\n",
      "lat = dataset.variables[\"lat\"][:]\n",
      "dataset.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "http://opendap.deltares.nl/thredds/dodsC/opendap/deltares/delwaq/NZBLOOM_2008_fPPtot_g_m-2_d-1.nc\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "# determine start time of calculations by determining unit\n",
      "\n",
      "# create indices array with m,n coordinate which fall with selection of\n",
      "# ices rectangles\n",
      "# HACK\n",
      "nodata = lon[0][0]\n",
      "assert nodata < -4e8, \"missing data expected in first cell\"\n",
      "\n",
      "indices = []\n",
      "# speed this up....\n",
      "# Match up the latitude longitudes with the polygon...\n",
      "# I don't understand why we're using a polygon\n",
      "for m in range(len(lon)):\n",
      "    for n in range(len(lon[m])):\n",
      "        if lon[m][n] != nodata and lat[m][n] != nodata:\n",
      "            pt=(lon[m][n],lat[m][n])\n",
      "            if apoly.contains(Point(pt)):\n",
      "                indices.append([m,n])\n",
      "indices = np.array(indices)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "# compute the mean of each polygon intersect\n",
      "pool = Pool(processes=4)\n",
      "listofrows = pool.map(mean_concentration, [(url, indices) for url in urls[:4]])\n",
      "allrows = []\n",
      "[allrows.extend(rows) for rows in listofrows]\n",
      "\n",
      "# create dataframe of vals and datetime index\n",
      "dfvals = pandas.DataFrame(allrows, columns=[\"time\", \"yearmonth\", \"concentration\"])\n",
      "\n",
      "grouped = dfvals.groupby('yearmonth').mean()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "AssertionError",
       "evalue": "Wrong number of items passed (3 vs 0)",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-7-6df312f0260f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;31m# create dataframe of vals and datetime index\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mdfvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpandas\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mallrows\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"time\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"yearmonth\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"concentration\"\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     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0mgrouped\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdfvals\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgroupby\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'yearmonth'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/fedorbaart/.virtualenvs/main/lib/python2.7/site-packages/pandas/core/frame.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, data, index, columns, dtype, copy)\u001b[0m\n\u001b[1;32m    421\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    422\u001b[0m                 mgr = self._init_ndarray(data, index, columns, dtype=dtype,\n\u001b[0;32m--> 423\u001b[0;31m                                          copy=copy)\n\u001b[0m\u001b[1;32m    424\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    425\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/fedorbaart/.virtualenvs/main/lib/python2.7/site-packages/pandas/core/frame.pyc\u001b[0m in \u001b[0;36m_init_ndarray\u001b[0;34m(self, values, index, columns, dtype, copy)\u001b[0m\n\u001b[1;32m    548\u001b[0m             \u001b[0mcolumns\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_ensure_index\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    549\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 550\u001b[0;31m         \u001b[0mblock\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmake_block\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    551\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mBlockManager\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mblock\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    552\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/fedorbaart/.virtualenvs/main/lib/python2.7/site-packages/pandas/core/internals.pyc\u001b[0m in \u001b[0;36mmake_block\u001b[0;34m(values, items, ref_items)\u001b[0m\n\u001b[1;32m    481\u001b[0m         \u001b[0mklass\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mObjectBlock\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    482\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 483\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mklass\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitems\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mref_items\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mndim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    484\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    485\u001b[0m \u001b[0;31m# TODO: flexible with index=None and/or items=None\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/fedorbaart/.virtualenvs/main/lib/python2.7/site-packages/pandas/core/internals.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, values, items, ref_items, ndim)\u001b[0m\n\u001b[1;32m     31\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     32\u001b[0m             raise AssertionError('Wrong number of items passed (%d vs %d)'\n\u001b[0;32m---> 33\u001b[0;31m                                  % (len(items), len(values)))\n\u001b[0m\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_ref_locs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mAssertionError\u001b[0m: Wrong number of items passed (3 vs 0)"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}