{
 "metadata": {
  "name": "",
  "signature": "sha256:71307476480b648887661e8aeea7766fb450154ce6f1d0aa555badf07fbd02a9"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import mmi\n",
      "import numpy as np\n",
      "import shapely.geometry\n",
      "import requests\n",
      "import zmq\n",
      "import zmq.eventloop\n",
      "import osgeo.osr\n",
      "import shapely.geometry\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 42
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ctx = zmq.Context()\n",
      "models = requests.get('http://localhost:22222/models').json()\n",
      "model = models[0]\n",
      "node = \"localhost\"\n",
      "\n",
      "pubport = model[\"ports\"][\"PUB\"]\n",
      "sub = ctx.socket(zmq.SUB)\n",
      "addr = \"tcp://{node:s}:{port:d}\".format(node=node, port=pubport)\n",
      "sub.connect(addr)\n",
      "sub.setsockopt_string(zmq.SUBSCRIBE, u'')\n",
      "\n",
      "pullport = model[\"ports\"][\"PULL\"]\n",
      "push = ctx.socket(zmq.PUSH)\n",
      "addr = \"tcp://{node:s}:{port:d}\".format(node=node, port=pullport)\n",
      "push.connect(addr)\n",
      "pushstream = zmq.eventloop.zmqstream.ZMQStream(push)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "names = {\"xk\", \"yk\", \"flowelemnode\"}\n",
      "for var in names:\n",
      "    mmi.send_array(pushstream, metadata={\"get_var\": var})\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "arrs = {}\n",
      "for name in names:\n",
      "    A, metadata = mmi.recv_array(sub)\n",
      "    arrs[metadata[\"name\"]] = A\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arrs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "{'flowelemnode': array([[  780,   853,   781,    -1,    -1,    -1],\n",
        "        [  781,   853, 13112,    -1,    -1,    -1],\n",
        "        [  781, 13112, 13110,    -1,    -1,    -1],\n",
        "        ..., \n",
        "        [28792, 28797, 28796, 28798,    -1,    -1],\n",
        "        [28793, 28798, 28796, 28794,    -1,    -1],\n",
        "        [28780, 28783, 28784, 28786, 28782, 28781]], dtype=int32),\n",
        " 'xk': array([ 547994.85513183,  547854.78975898,  547712.6433233 , ...,\n",
        "         595202.07393365,  595315.15648816,  595225.4700929 ]),\n",
        " 'yk': array([ 4204745.95785394,  4204848.64553312,  4204954.42679592, ...,\n",
        "         4214406.36377765,  4214530.90081606,  4214617.57328289])}"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flowelemnode = np.ma.masked_equal(arrs[\"flowelemnode\"], -1)\n",
      "xk = arrs[\"xk\"]\n",
      "yk = arrs[\"yk\"]\n",
      "\n",
      "points = np.c_[xk, yk]\n",
      "src_srs = osgeo.osr.SpatialReference()\n",
      "src_srs.ImportFromEPSG(model[\"epsg\"])\n",
      "dst_srs = osgeo.osr.SpatialReference()\n",
      "dst_srs.ImportFromEPSG(4326)\n",
      "\n",
      "\n",
      "\n",
      "transform = osgeo.osr.CoordinateTransformation(src_srs, dst_srs)\n",
      "wkt_points = np.array(transform.TransformPoints(points))\n",
      "wkt_points"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 47,
       "text": [
        "array([[-122.45342788,   37.98909466,    0.        ],\n",
        "       [-122.45501605,   37.99002753,    0.        ],\n",
        "       [-122.45662776,   37.99098838,    0.        ],\n",
        "       ..., \n",
        "       [-121.91460836,   38.0724204 ,    0.        ],\n",
        "       [-121.91330272,   38.07353071,    0.        ],\n",
        "       [-121.91431353,   38.07432118,    0.        ]])"
       ]
      }
     ],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# compute number of nodes per element\n",
      "elemtype = (~flowelemnode.mask).sum(1)\n",
      "# compute the places where we can split the 1d array\n",
      "splitidx = np.cumsum(np.r_[elemtype][:-1])\n",
      "# flatten the node numbers, remove the missings\n",
      "idx = flowelemnode[(~flowelemnode.mask)]-1\n",
      "# use the indices to lookup lon, lat\n",
      "# only keep x and y, transform to list\n",
      "polycoords = [x[:,:2] for x in np.split(wkt_points[idx],splitidx)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "polys = [shapely.geometry.mapping(shapely.geometry.Polygon(xy)) for xy in polycoords]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 70
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "features = [dict(type=\"Feature\", geometry=poly, properties=dict(index=i)) for i, poly in enumerate(polys)]\n",
      "collection = dict(type=\"FeatureCollection\", features=features)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 71
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import json\n",
      "json.dump(collection, open(\"grid.json\",\"w\"))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 72
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 68,
       "text": [
        "[array([[-122.4806033 ,   38.09635956],\n",
        "        [-122.47300867,   38.10211462],\n",
        "        [-122.48402421,   38.09970107]]),\n",
        " array([[-122.48402421,   38.09970107],\n",
        "        [-122.47300867,   38.10211462],\n",
        "        [-122.48052373,   38.10561733]]),\n",
        " array([[-122.48402421,   38.09970107],\n",
        "        [-122.48052373,   38.10561733],\n",
        "        [-122.48459261,   38.10532472]]),\n",
        " array([[-122.48402421,   38.09970107],\n",
        "        [-122.48459261,   38.10532472],\n",
        "        [-122.48897247,   38.1033987 ]]),\n",
        " array([[-122.48897247,   38.1033987 ],\n",
        "        [-122.48459261,   38.10532472],\n",
        "        [-122.48771973,   38.10716098]]),\n",
        " array([[-122.47300867,   38.10211462],\n",
        "        [-122.46462764,   38.1077655 ],\n",
        "        [-122.4757648 ,   38.10873706]]),\n",
        " array([[-122.47300867,   38.10211462],\n",
        "        [-122.4757648 ,   38.10873706],\n",
        "        [-122.48052373,   38.10561733]]),\n",
        " array([[-122.46462764,   38.1077655 ],\n",
        "        [-122.47193063,   38.11491708],\n",
        "        [-122.4757648 ,   38.10873706]]),\n",
        " array([[-122.46462764,   38.1077655 ],\n",
        "        [-122.45617388,   38.11251683],\n",
        "        [-122.45986205,   38.11667626]]),\n",
        " array([[-122.46462764,   38.1077655 ],\n",
        "        [-122.45986205,   38.11667626],\n",
        "        [-122.47193063,   38.11491708]])]"
       ]
      }
     ],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}