Friday, November 16, 2012

Projecting some thoughts about pyproj.

So I've been playing around with pyproj a little lately and have come up with a couple new thoughts.  As always I'm typing as I go so not everything is a clean example of good code or procedure, just something I've run as a test of my thoughts.

First, the application; below is which when run displays the lat long of the current cursor position in the DspMsg0 area (the top long line right below the command area).  I decided to do this as an app so it is running in real time and clicking the buttons just performs specific tasks.  Speaking of the buttons, B1 or the left mouse button prints the current position in the local coordinates, long lat, and notes the projection used. Something like

1369447.29077 471184.7857 (-84.66858719560254, 40.93954025528449) epsg:3753

In this part of the world we primarily deal with two projections so I kept it simple and just used B2 or the right mouse button to toggle between the two.  It would have been easy enough to bring up some kind of selection dialog, but for now I only need the two. The current projection is also displayed in the DspMsg area.

I recently discovered why I always had to send and receive meters even though the EPSG code clearly specified at foot based projection.  It was by design so there was never any question what was expected or returned.  The answer was always meters, however the code has since been modified to include an init keyword


which handily enough, does exactly what it says. I understand the value of the consistency, but appreciate the option. I have also discovered in testing that you can take a proj4 definition and paste it straight into an object init function in quotes and get the exact same functionality as using the EPSG codes.  For example...

ohs=pyproj.Proj('+proj=lcc +lat_1=40.03333333333333 +lat_2=38.73333333333333 +lat_0=38 +lon_0=-82.5 +x_0=600000 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs' )

Well in any case here is the app, it isn't worth much except by way of example, but pyproj itself is really handy when working with things like flight prep or control layouts when you want to give somebody geodetic coordinates to work from.

print 'showll modified 12:17 PM 11/14/2012'
# Show Longitude Latitude of current position.
Copyright 2012 Dennis Shimer
Vr Mapping copyright Cardinal Systems, LLC
No warranties as to safety or usability expressed or implied.
Free to use, copy, modify, distribute with author credit.

Only two projections are supported. Both are NAD83(HARN) / Ohio (ftUS).
EPSG:3753 for North Zone
EPSG:3754 for South Zone
Button 2 (right mouse) toggles the projection, Button 1 (left mouse)
will print out the current position in local coordinate and it
Longitude, Latitude decimal degrees.
MYAPPID = 10011

Gui = PyVrGui()

import pyproj,sys

global Projection

Projection=pyproj.Proj(init='EPSG:3754', preserve_units=True)
Gui.DspMsg('OH S EPSG:3754')

def MyAppDigCB (x, y, z, key, id):

 global Projection
 if (key == 1):
  print x,y,repr(Projection(x,y,inverse=True)),Projection.srs.split('=')[1]
 elif (key == 2):
  if Projection.srs.count('3754'):
   Projection=pyproj.Proj(init='EPSG:3753', preserve_units=True)
   Gui.DspMsg('OH N EPSG:3753')
  elif Projection.srs.count('3753'):
   Projection=pyproj.Proj(init='EPSG:3754', preserve_units=True)
   Gui.DspMsg('OH S EPSG:3754')
 elif (key == 11):
  MyApp.Close ()
  MyApp.StopRunning ()

def MyAppKeyCB (key):

 MyAppDigCB (x,y,z, key, 0)

 MyApp = PyVrApp (MYAPPID, "ShowLL")
 print "Could not start application:", sys.exc_info()[0]
 Mk = PyVrMenuKeys ("ShowLL", 10)
 Mk.SetMkLabels ("ShowLL","1 Print pt","2 Tog N/S", "3", "4", "5", "6", "7", "8", "9", "*", "0", "# Quit")
 print "Could not create menukeys:", sys.exc_info()[0]

 MyApp.SetDigCB (MyAppDigCB)
 Mk.SetKeyCB (MyAppKeyCB)
 MyApp.KeepRunning ()
 Mk.Delete ()
 print "Exception during application:", sys.exc_info()[0]
 MyApp.Close ()
 MyApp.StopRunning ()

How about one more bonus snippet with regard just switching from on projection to another. Thanks to Proj4 and pyproj, this is just how easy it can be.

>>> ohs=pyproj.Proj(init='EPSG:3754', preserve_units=True)
>>> ohn=pyproj.Proj(init='EPSG:3753', preserve_units=True)
>>> pyproj.transform(ohs,ohn,1759595.61,661026.83)
(1759551.3008025475, 54039.04746006191)

No comments:

For anyone interested in trying VrPython for the first time or if you are early in the game, I suggest going to the earliest posts and working forward. I use VrPython every day for many wonderful things, needless to say it will change and could potentially damage a file. Any risk associated with using VrPython or any code or scripts mentioned here lies solely with the end user.

The "Personal VrPython page" in the link section will contain many code examples and an organized table of contents to this blog in a fairly un-attractive (for now) form.