{
 "metadata": {
  "name": "",
  "signature": "sha256:3cb8e8201e238e76c39b6f53a5379a8e077dc0d971d25f6903d35aab15d938d8"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# how to convert Dutch coordinate system epsg:28992 with order 20 cm accuracy\n",
      "\n",
      "# https://rdinfo.kadaster.nl/pics/publijst1.gif\n",
      "# https://rdinfo.kadaster.nl/?inhoud=/rd/info.html%23publicatie&navig=/rd/nav_serverside.html%3Fscript%3D1\n",
      "import haversine # for differences\n",
      "D = {}\n",
      "D[\"Puntnummer\"]        = '019111'\n",
      "D[\"Actualiteitsdatum\"] = '1999,6,1'\n",
      "D[\"Nr\"]                = 17\n",
      "D[\"X\"]                 = 155897.26\n",
      "D[\"Y\"]                 = 603783.39\n",
      "D[\"H\"]                 = 3.7\n",
      "D[\"lat\"]               = 53+25/60+13.2124/3600\n",
      "D[\"lon\"]               =  5+24/60+02.5391/3600\n",
      "D[\"h\"]                 = 44.83"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# USING EPSG\n",
      "# Dutch RD: this DID not work in the past not work due to an error in proj4 and the epsg database\n",
      "# (http://trac.osgeo.org/geotiff/ticket/22)\n",
      "\n",
      "srs1 = pyproj.Proj(init='epsg:4326' ) # wgs84\n",
      "srs2 = pyproj.Proj(init='epsg:28992')\n",
      "x,y  = pyproj.transform(srs1, srs2, D[\"lon\"]+0.,D[\"lat\"]+0.) # Do add 0. to avoid trunc issues. \n",
      "print('x:',D[\"X\"],x,D[\"X\"]-x)\n",
      "print('y:',D[\"Y\"],y,D[\"Y\"]-y)\n",
      "print('error ',sqrt((D[\"Y\"]-y)**2 +(D[\"X\"]-x)**2),'m')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "x: 155897.26 155897.35463811128 -0.09463811127352528\n",
        "y: 603783.39 603783.5452836686 -0.15528366854414344\n",
        "error  0.181849910151 m\n"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# USING PROJ4 string\n",
      "import pyproj\n",
      "proj4_params = '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +no_defs'\n",
      "\n",
      "srs1 = pyproj.Proj(init='epsg:4326' ) # wgs84\n",
      "srs2 = pyproj.Proj(proj4_params)\n",
      "x,y  = pyproj.transform(srs1, srs2, D[\"lon\"]+0.,D[\"lat\"]+0.) # Do add 0. to avoid trunc issues. \n",
      "print('x:',D[\"X\"],x,D[\"X\"]-x)\n",
      "print('y:',D[\"Y\"],y,D[\"Y\"]-y)\n",
      "print('error ',sqrt((D[\"Y\"]-y)**2 +(D[\"X\"]-x)**2),'m')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "x: 155897.26 155897.35478599914 -0.09478599912836216\n",
        "y: 603783.39 603783.5583675711 -0.16836757108103484\n",
        "error  0.193214969975 m\n"
       ]
      }
     ],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#import osgeo\n",
      "#s_srs  = osgeo.osr.SpatialReference()\n",
      "#s_srs.ImportFromEPSG(28992 )\n",
      "#lon,lat = sp_srs(X+0.,Y+0.,inverse=True) # inverse means x/y to lat/lon. Thsi does not work due to error in EPSG database. Do add 0. to avoid trunc issues\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}