{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Let's check our version\n",
      "!cs2cs \n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rel. 4.7.1, 23 September 2009\r\n",
        "usage: cs2cs [ -eEfIlrstvwW [args] ] [ +opts[=arg] ]\r\n",
        "                   [+to [+opts[=arg] [ files ]\r\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# We would expect this to work, but it doesn't\n",
      "# The towgs84 is used for transformation\n",
      "! echo \"49554 215020 0\" | cs2cs -v +init=epsg:31370 +to +init=epsg:4326  -f %.6f"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "# ---- From Coordinate System ----\r\n",
        "#Lambert Conformal Conic\r\n",
        "#\tConic, Sph&Ell\r\n",
        "#\tlat_1= and lat_2= or lat_0\r\n",
        "# +init=epsg:31370 +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339\r\n",
        "# +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438\r\n",
        "# +ellps=intl\r\n",
        "# +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m\r\n",
        "# +no_defs\r\n",
        "# ---- To Coordinate System ----\r\n",
        "#Lat/long (Geodetic alias)\r\n",
        "#\t\r\n",
        "# +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\r\n",
        "# +towgs84=0,0,0\r\n",
        "#--- following specified but NOT used\r\n",
        "# +ellps=WGS84\r\n",
        "2.928019\t51.235804 347.868509\r\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# We get a different result if we use invproj, what's going on?\n",
      "! echo \"49554 215020\" | invproj -v +init=epsg:31370 -f %.6f\n",
      "\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "#Lambert Conformal Conic\r\n",
        "#\tConic, Sph&Ell\r\n",
        "#\tlat_1= and lat_2= or lat_0\r\n",
        "# +init=epsg:31370 +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339\r\n",
        "# +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438\r\n",
        "# +ellps=intl\r\n",
        "# +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m\r\n",
        "# +no_defs\r\n",
        "2.929249\t51.236883\r\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# We get the correct result if we fix the last parameter of the towgs84 parameter\n",
      "# http://osgeo-org.1560.x6.nabble.com/proj4-epsg-31370-possible-wrong-parameters-td3841364.html\n",
      "!echo \"49554 215020\" | cs2cs -v \\\n",
      "    +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 \\\n",
      "    +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl \\\n",
      "    +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747 +units=m \\\n",
      "    +no_defs -to  +proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs -f %6f"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "o ---- From Coordinate System ----\r\n",
        "#Lambert Conformal Conic\r\n",
        "#\tConic, Sph&Ell\r\n",
        "#\tlat_1= and lat_2= or lat_0\r\n",
        "# +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90\r\n",
        "# +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl\r\n",
        "# +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747\r\n",
        "# +units=m +no_defs +datum=WGS84\r\n",
        "#--- following specified but NOT used\r\n",
        "# +proj=longlat +ellps=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0\r\n",
        "o ---- To Coordinate System ----\r\n",
        "#Lat/long (Geodetic alias)\r\n",
        "#\t\r\n",
        "# +proj=latlong +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0\r\n",
        "2.930480\t51.236358 41.428643\r\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# This does not work for invproj because proj and invproj ignore the datum shift\n",
      "# The towgs84 is ignored\n",
      "!echo \"49554 215020\" | invproj -v \\\n",
      " +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 \\\n",
      " +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl \\\n",
      "    +towgs84=-1,1,1,1,1,1,1 +units=m \\\n",
      " +no_defs \\\n",
      " -f %6f"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "#Lambert Conformal Conic\r\n",
        "#\tConic, Sph&Ell\r\n",
        "#\tlat_1= and lat_2= or lat_0\r\n",
        "# +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90\r\n",
        "# +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl\r\n",
        "# +towgs84=-1,1,1,1,1,1,1 +units=m +no_defs\r\n",
        "2.929249\t51.236883\r\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# There is another datums shift that gives the same result on spatialreference.org\n",
      "# http://spatialreference.org/ref/sr-org/epsg31370-correct-datum-shift/\n",
      "!echo \"49554 215020\" | cs2cs -v \\\n",
      "+proj=lcc +lat_1=51.16666723333334 +lat_2=49.83333389999999 \\\n",
      "+lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 \\\n",
      "+ellps=intl +towgs84=-99.1,53.3,-112.5,0.419,-0.83,1.885,-1.0 +units=m +no_defs \\\n",
      "-to +datum=WGS84 +ellps=WGS84 -f %6f"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "o ---- From Coordinate System ----\r\n",
        "#Lambert Conformal Conic\r\n",
        "#\tConic, Sph&Ell\r\n",
        "#\tlat_1= and lat_2= or lat_0\r\n",
        "# +proj=lcc +lat_1=51.16666723333334 +lat_2=49.83333389999999 +lat_0=90\r\n",
        "# +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl\r\n",
        "# +towgs84=-99.1,53.3,-112.5,0.419,-0.83,1.885,-1.0 +units=m +no_defs\r\n",
        "# +datum=WGS84\r\n",
        "#--- following specified but NOT used\r\n",
        "# +ellps=WGS84 +ellps=WGS84 +towgs84=0,0,0\r\n",
        "o ---- To Coordinate System ----\r\n",
        "#Lat/long (Geodetic alias)\r\n",
        "#\t\r\n",
        "# +proj=latlong +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0\r\n",
        "2.930478\t51.236358 41.262401\r\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can also use the ogr2ogr tool, which uses the osr library.\n",
      "\n",
      "This one is more up to date on my system. We have to define the format of our csv file in a vrt file (or use geojson).\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file test.vrt\n",
      "<OGRVRTDataSource>\n",
      "    <OGRVRTLayer name=\"test\">\n",
      "        <SrcDataSource>test.csv</SrcDataSource>\n",
      "        <GeometryType>wkbPoint</GeometryType>\n",
      "        <LayerSRS>EPSG:31370</LayerSRS>\n",
      "        <GeometryField encoding=\"PointFromColumns\" x=\"x\" y=\"y\"/>\n",
      "    </OGRVRTLayer>\n",
      "</OGRVRTDataSource>"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Overwriting test.vrt\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "! echo \"x,y\\n49554,215020\"  >test.csv \n",
      "jsontxt = ! ogr2ogr -f GeoJSON -t_srs EPSG:4326 /dev/stdout test.vrt\n",
      "import json\n",
      "#json.loads(jsontxt)\n",
      "print(\"\\n\".join(jsontxt))\n",
      "json.loads(\"\".join(jsontxt))['features'][0]['geometry']['coordinates']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "{\n",
        "\"type\": \"FeatureCollection\",\n",
        "\"features\": [\n",
        "{ \"type\": \"Feature\", \"properties\": { \"x\": \"49554\", \"y\": \"215020\" }, \"geometry\": { \"type\": \"Point\", \"coordinates\": [ 2.930479615690256, 51.236357952698896 ] } }\n",
        "\n",
        "]\n",
        "}\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "[2.930479615690256, 51.236357952698896]"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Or you can use a python script to talk to the osr library, which as the correct coordinates.\n",
      "import osgeo.osr\n",
      "src = osgeo.osr.SpatialReference()\n",
      "dst = osgeo.osr.SpatialReference()\n",
      "src.ImportFromEPSG(31370)\n",
      "print src.GetTOWGS84()\n",
      "dst.ImportFromEPSG(4326)\n",
      "transform = osgeo.osr.CoordinateTransformation(src, dst)\n",
      "lon, lat, z = transform.TransformPoint(49554, 215020, 0)\n",
      "print(lon, lat, z)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(-106.869, 52.2978, -103.724, 0.3366, -0.457, 1.8422, -1.2747)\n",
        "(2.930479615690256, 51.236357952698896, 41.428643393330276)\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    }
   ],
   "metadata": {}
  }
 ]
}