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)

Tuesday, October 23, 2012

Where in the world... (and by that I mean geographically, not disappearing for 3 years)

Honestly,  first post in 3 years?  Well yes in fact it is, not much to say and probably nobody that really cares so it is what it is.  Been busy coding and there are some truly awesome libraries out there.  Lots of work extracting information out of Vr and dropping it directly into XLS spreadsheets using xlwt and of course the switch over to 2.7 (and 64 bit) has provided plenty of interesting experience. However, what is so important that I'm actually going to sit and type for a few minutes? It is pyproj of course.  Up to now I had done all my re-projection using the corpscon dll libraries which to date I can not find any way of porting to 2.7 64bit.  In the mean time I have learned a lot more about great open libraries like PROJ.4 but just never had the time to investigate.  Necessity being the mother of invention (or in this case studying) I needed a replacement for CorpsCon and the answer was pyproj.  A snippet is worth a million words so here goes.

After downloading and installing pyproj I start by importing the library
import pyproj

Next we create a Projection object and initialize it with a particular projection in this case Ohio State Plane North 1938 EPSG code 2834 

If I have an X,Y coordinate I can just allow the Projection object to compute the Long, Lat.  Notice that the default must be to go from geodetic to X,Y so the "inverse=True" is there and also that it is expecting (and returns) meters so we need to convert.  Both these statements may be a little off but I've only worked with the library for about 30 minutes so far and haven't actually read the docs.
(-83.84694414368035, 41.58018579703697)

So what we get back is a list with X and Y.  If I had passed in a Long and Lat position it would again return a list of coordinates in meters so lets create a quick function to take the coordinate list and return feet.
def to_ft((x,y)):
    return x*3.280833333333333,y*3.280833333333333

Then simply do the projection without the inverse.   
(1830492.2997760142, 486159.352289076)

In the docs there are easy ways to create two different Projection objects and go directly from one to the other, but I pretty much just needed to go from state plane to geodetic as the basis for creating KML output. The change in libraries dropped about 75 lines of code, you gotta love python libraries.

Speaking of libraries, I found a site that hosts every python library I could think of (and about a hundred I couldn't) as a binary Windows installer in everything from Python 2.5 32bit to Python 3.2 64bit.  Find them at

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.