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 
Projection=pyproj.Proj(init='EPSG: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.
Projection(1600000/3.280833333,700000/3.280833333,inverse=True)
(-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.   
to_ft(Projection(-83.0,41.0))
(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 http://www.lfd.uci.edu/~gohlke/pythonlibs/

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.