convert_wgs83_to_utm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# convert decimal lat / lon in WGS83 EPSG: 4326 to NAD83(NRS2007) UTM meters, Zone 17N EPSG:3724\n",
    "# to get this to work, I had to define GDAL_DATA in my Win environment path\n",
    "\n",
    "# csherwood@usgs.gov\n",
    "# 4 November 2016"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Disclaimer\n",
    "\n",
    "Unless otherwise noted, the software and content on this site is\n",
    "in the public domain because it contains materials developed by\n",
    "the United States Geological Survey, an agency of the United States\n",
    "Department of Interior. For more information, see the official USGS\n",
    "copyright policy at:\n",
    "\n",
    "http://www.usgs.gov/visual-id/credit_usgs.html#copyright\n",
    "\n",
    "This software and content is distributed on an \"AS IS\" BASIS, WITHOUT\n",
    "WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. Any\n",
    "dependent libraries included here are distributed under open source\n",
    "(or open source-like) licenses/agreements. Appropriate license agreements\n",
    "are included with each library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " >> SOURCE COORDINATE REFERENCE SYSTEM: \n",
      "GEOGCS[\"WGS 84\",\n",
      "    DATUM[\"WGS_1984\",\n",
      "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
      "            AUTHORITY[\"EPSG\",\"7030\"]],\n",
      "        AUTHORITY[\"EPSG\",\"6326\"]],\n",
      "    PRIMEM[\"Greenwich\",0,\n",
      "        AUTHORITY[\"EPSG\",\"8901\"]],\n",
      "    UNIT[\"degree\",0.0174532925199433,\n",
      "        AUTHORITY[\"EPSG\",\"9122\"]],\n",
      "    AUTHORITY[\"EPSG\",\"4326\"]]\n",
      " >> DESTINATION COORDINATE REFERENCE SYSTEM: \n",
      "PROJCS[\"NAD83(NSRS2007) / UTM zone 17N\",\n",
      "    GEOGCS[\"NAD83(NSRS2007)\",\n",
      "        DATUM[\"NAD83_National_Spatial_Reference_System_2007\",\n",
      "            SPHEROID[\"GRS 1980\",6378137,298.257222101,\n",
      "                AUTHORITY[\"EPSG\",\"7019\"]],\n",
      "            TOWGS84[0,0,0,0,0,0,0],\n",
      "            AUTHORITY[\"EPSG\",\"6759\"]],\n",
      "        PRIMEM[\"Greenwich\",0,\n",
      "            AUTHORITY[\"EPSG\",\"8901\"]],\n",
      "        UNIT[\"degree\",0.0174532925199433,\n",
      "            AUTHORITY[\"EPSG\",\"9122\"]],\n",
      "        AUTHORITY[\"EPSG\",\"4759\"]],\n",
      "    PROJECTION[\"Transverse_Mercator\"],\n",
      "    PARAMETER[\"latitude_of_origin\",0],\n",
      "    PARAMETER[\"central_meridian\",-81],\n",
      "    PARAMETER[\"scale_factor\",0.9996],\n",
      "    PARAMETER[\"false_easting\",500000],\n",
      "    PARAMETER[\"false_northing\",0],\n",
      "    UNIT[\"metre\",1,\n",
      "        AUTHORITY[\"EPSG\",\"9001\"]],\n",
      "    AXIS[\"Easting\",EAST],\n",
      "    AXIS[\"Northing\",NORTH],\n",
      "    AUTHORITY[\"EPSG\",\"3724\"]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# GDAL Python library\n",
    "from osgeo import osr\n",
    "# EXIF Reader\n",
    "import exifread\n",
    "# Image library\n",
    "#from PIL import Image\n",
    "\n",
    "from IPython.display import Image\n",
    "Image(\"vdatum_screen_grab.jpg\")\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "get_ipython().magic(u'matplotlib inline')\n",
    "\n",
    "src_crs = 'EPSG:4326'\n",
    "dst_crs = 'EPSG:3724'\n",
    "\n",
    "# source coordinate system\n",
    "srs_cs = osr.SpatialReference()\n",
    "srs_cs.SetFromUserInput(src_crs)\n",
    "\n",
    "# destination coordinate system\n",
    "dst_cs = osr.SpatialReference()\n",
    "dst_cs.SetFromUserInput(dst_crs)\n",
    "\n",
    "# osr image transformation object\n",
    "transform = osr.CoordinateTransformation(srs_cs, dst_cs)\n",
    "\n",
    "# print coordinate system information\n",
    "print \" >> SOURCE COORDINATE REFERENCE SYSTEM: \"\n",
    "print srs_cs\n",
    "\n",
    "print \" >> DESTINATION COORDINATE REFERENCE SYSTEM: \"\n",
    "print dst_cs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "    # I checked this transform using VDATUM for the Point 1 south below, with source:\n",
    "    # Source: \"WGS84(transit) - use NAD83 (see NGS's HTDP)\"\n",
    "    # Target: \"NAVD88\" using GEOID12B\n",
    "    # northing and easting transforms were reproduced exactly\n",
    "    # for Point 1 south, height was 1.0326 m NADV88, with vertical uncertainty of 7.6158 cm\n",
    "    # deltaz = 28.7909\n",
    "    # for Point 9 north, deltaz - 28.7911\n",
    "    # so, using:\n",
    "    deltaz = 28.7910"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Label,Lat (degrees),Lat (minutes),Lat (seconds),Lon (degrees),Lon (minutes),Lon (seconds),Z (Height above the WGS84 ellipsoide),Lat (WGS84),Lon(WGS84)\n",
      "\n",
      "Point 1 south, 29.683666, -81.218714, -27.7583, 478838.9754, 3283753.8587, 1.0327\n",
      "Point2 south, 29.683818, -81.218888, -28.6757, 478822.1098, 3283770.7518, 0.1153\n",
      "Point3 south, 29.683666, -81.218937, -28.7528, 478817.3594, 3283753.8340, 0.0382\n",
      "Point4 south, 29.683769, -81.218988, -28.5370, 478812.4735, 3283765.2579, 0.2540\n",
      "Point5 south, 29.683624, -81.218852, -28.1598, 478825.6177, 3283749.2133, 0.6312\n",
      "Point6 south, 29.683519, -81.218881, -27.6961, 478822.7810, 3283737.5846, 1.0949\n",
      "Point7 south, 29.683355, -81.218757, -27.5410, 478834.6762, 3283719.3903, 1.2500\n",
      "Point8 south, 29.683599, -81.218536, -28.6292, 478856.1846, 3283746.3694, 0.1618\n",
      "Point 1 north, 29.685358, -81.219510, -27.7330, 478762.2574, 3283941.4009, 1.0580\n",
      "Point 2 north, 29.685358, -81.219427, -27.5740, 478770.2724, 3283941.4085, 1.2170\n",
      "Point 3 north, 29.685595, -81.219426, -27.9800, 478770.4969, 3283967.6541, 0.8110\n",
      "Point 4 north, 29.685327, -81.219665, -28.9210, 478747.2610, 3283938.0370, -0.1300\n",
      "Point 5 north, 29.685677, -81.219520, -28.1661, 478761.3514, 3283976.7921, 0.6249\n",
      "Point 6 north, 29.686064, -81.219693, -28.0038, 478744.7407, 3284019.7250, 0.7872\n",
      "Point 7 north, 29.686207, -81.219834, -27.9100, 478731.1534, 3284035.5378, 0.8810\n",
      "Point 8 north, 29.686231, -81.219778, -27.7394, 478736.5016, 3284038.2464, 1.0516\n",
      "Point 9 north, 29.686201, -81.219734, -28.0870, 478740.8269, 3284034.8617, 0.7040\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<function close>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fnin = \"D:\\\\crs\\\\proj\\\\2017_Matthew\\\\maitane\\\\2016-10-30_control_points.csv\"\n",
    "fnout = \"D:\\\\crs\\\\proj\\\\2017_Matthew\\\\maitane\\\\2016-10-30_GCP_UTM.csv\"\n",
    "\n",
    "header_out = \"Label, Lat (WGS84), Lon (WGS84), Alt. (m), Easting (m NAD83 UTM17N), Northing (m NAD83 UTM17N), Elev. (m NADV88)\"\n",
    "f=open(fnin,'r')\n",
    "fo = open(fnout, 'w')\n",
    "header = f.readline()\n",
    "print header\n",
    "fo.write(header_out+\"\\n\")\n",
    "for line in f:\n",
    "    line = line.strip()\n",
    "    col = line.split(',')\n",
    "    name = col[0]\n",
    "    lat = float(col[1])+float(col[2])/60.+float(col[3])/3600.\n",
    "    lon = -(float(col[4])+float(col[5])/60.+float(col[6])/3600.)\n",
    "    ze = float(col[7])\n",
    "    e,n,zt = transform.TransformPoint(lon, lat, ze)\n",
    "    znavd88 = zt+deltaz\n",
    "    #print(name, lat, lon, ze, e, n, znavd88)\n",
    "    print(\"{0:s}, {1:.6f}, {2:.6f}, {3:.4f}, {4:.4f}, {5:.4f}, {6:.4f}\".format(name,lat,lon,ze,e,n,znavd88))\n",
    "    fo.write(\"{0:s}, {1:.6f}, {2:.6f}, {3:.4f}, {4:.4f}, {5:.4f}, {6:.4f}\\n\".format(name,lat,lon,ze,e,n,znavd88))\n",
    "f.close\n",
    "fo.close"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}