Almost sometimes biennial proof of life.

 Recent conversation about a particular legend in a particular field went something like "Did you realize he was that old?" and my first thought was "I didn't even realize he was still alive!".  Makes me wonder how many people think that when the look at such a long-ignored blog.  Well I'm not that old and I am still alive and the question I was asked was "Are you still writing any python?", to which the answer is yes almost daily but I also realize that no one really cares.

I did pass off my "BeginnersVrPython3.py" script to someone and realize that I should make a quick video to document what it does in case anyone in the office wanted to start plinking on coding some python.

It can be found at

https://rumble.com/v6ezy07-cardinal-systems-vrone-vrtwo-08-beginning-vrpython.html

and the script can be found at

https://drive.google.com/file/d/1Kyt0bkb0qPEAgtl7HriiFfjVmt7s0dzh/view?usp=sharing

Feel free to contact me at the email in the script if you have any questions.

Placing Lines Within Polygons Uniquely

 I'll be honest, there might be a function to do this but I couldn't find it and it was only going to take a few minutes so let's take a look.  I have a bunch of triangles for and  XML file that need to be subdivided into areas that are defined by polygons.  I need to make sure that the triangles are kept intact and that there are no duplicates.  In my case it didn't need to be an exact science and I wasn't worried about the edges.  My solution was to check each triangle against a polygon and if the first point of the triangle was within the polygon it would be flagged, or in this case the layer changed.  I just banged it out and I may make changes but for now it does the trick.

print ('line_change_line_in_polygon.py modified   06:25 2023/11/07')

# Change layer of line that starts within polygon.

'''

Author: Dennis Shimer dshimer@gmail.com

License: cc-by-sa  http://creativecommons.org/licenses/by-sa/3.0/

'''

Ws=PyVrWs()

Line=PyVrLine()

BoundingLine=PyVrLine()

Gui=PyVrGui()

CheckLayer=600

NewLayer=1001

PointsFound=0

WsNum=Ws.Aws()

PromBox = VrPromBox ('Check lines in Polygon.', 30, 1)

PromBox.AddInt ('Layer to check', CheckLayer, 1, 32000)

PromBox.AddInt ('Change layer to', NewLayer, 1, 32000)

if (PromBox.Display(0) == 0):

CheckLayer=PromBox.GetIntByPrompt('Layer to check')

NewLayer=PromBox.GetIntByPrompt('Change layer to')

BoundingLine.Id()

print (CheckLayer,BoundingLine.GetLayer())

Gui.ProgInit('Lines',Ws.GetLineCount(WsNum))

for EntNum in range ( Ws.GetLineCount(WsNum)):

Gui.ProgSet(EntNum)

if CheckLayer == Ws.GetLineLayer(WsNum,EntNum):

Line.Load (WsNum, EntNum)

X,Y,Z=Line.GetPoint(0)

if BoundingLine.IsPointInside (X,Y):

# print (X,Y,Z)

PointsFound=PointsFound+1

Line.SetLayer(NewLayer)

Line.ReRec(1)

# Line.Plot()

Gui.ProgReset()

print (PointsFound)

 

Modifying LiDAR Intensity in Vr Applications

I actually don't have the patience to look back through posts with regard to modifying LiDAR data so I'll just drop this one in as an example of how simple it can be.  This will have limited usefulness because it is fixing an apparent shortcoming in Vr that will no doubt be rectified at any moment, but hopefully the example will continue to be useful.  The current problem is that points created via dsmare (DSM Area) have a nice RGB value, but sometimes it would be nice to have Intensity as well. Rather than making it tricky I found that if I just add the R, G, and B values and divide by 3 it provides a useful "intensity like" display.    
There are only a couple lines that do anything productive so I'll just explain them.

print 'rgb2int.py modified 12:23 PM 9/25/2018'
# Set point intensity based on rgb
'''
License: cc-by-sa  http://creativecommons.org/licenses/by-sa/3.0/

'''
Ws=PyVrWs()
Punt=PyVrPunt()
Gui=PyVrGui()
WsNum=Ws.Aws()
PointBufferCount=Ws.GetPuntBufCount(WsNum)

Gui.ProgInit('Points',Ws.GetPuntBufCount(WsNum))
for PointBufferNumber in range(PointBufferCount):
    Punt.Load(WsNum,PointBufferNumber)
    Gui.ProgSet(PointBufferNumber)
    for PuntNum in range(Punt.GetCount()):
        x,y,z,PuntA,r,g,b= Punt.GetPuntRgb (PuntNum)
        PuntA['Int']=int((r+g+b)/3)
        Punt.ChgPuntA (PuntNum, PuntA)
    Punt.ReRec()
Gui.ProgReset()
PyVrGr().Replot()

 You need to have a PyVrPunt object to work with.  This object will be used to load up a "Point Buffer" full of LiDAR data which we can then step through the individual point data.  Note that you need to loop over the whole file loading each of the point buffers into the Punt object. This is done with the first for loop

for PointBufferNumber in range(PointBufferCount):

As we step through all the buffers, each one needs to be loaded and then it will be stepped through loading up the individual points hence the next for loop.

for PuntNum in range(Punt.GetCount()):

There are several ways to get particular data items from individual points and much of the data is stored in an attribute dictionary known as PuntA.  Many times I just load up PuntA, query it, make decisions, then make changes based on what I find.  In this case it is handy that the method that extracts the R,G,B data also pulls out PuntA at the same time because Intensity is stored within PuntA.  The first thing that needs doing is grabbing the data for the individual point.


x,y,z,PuntA,r,g,b= Punt.GetPuntRgb (PuntNum)

Then we just turn right around and push the intensity into the current PuntA dictionary and record the modified attributes back to the point.


PuntA['Int']=int((r+g+b)/3)
Punt.ChgPuntA (PuntNum, PuntA)

Once all the points are changed in the currently loaded buffer (Punt), the whole buffer can be re-recorded.


Punt.ReRec()

The Punt class is fully documented in the Python Programming section of the documentation but it basically comes down to knowing if you are creating or modifying something inside the attribute dictionary or outside it and working with the appropriate methods.

Musings on being a senior programmer

In some context's the title would mean being an experienced, admired, authoritative person.  Alas for me it just means being old, and in retrospect using the word programmer is a bit of a stretch.  Maybe it is harder to remember things because every day I have so much more important information to pack away for later retrieval, or maybe it is just harder to remember.  As we evolve as a company I find myself adding layer names and function keys that I know I'll use again but not often enough that I'll probably remember them.  Gone are the days when a function was one of 120 three digit numbers (which of course I still remember because I've been using most of them for 30 years). I quit trying to remember things I can Google long ago so why not treat the new layers and functions the same, after all if I expose myself to them often enough I'll eventually remember, and if I don't then searching is easy enough.  Enter sealay and seafk.  Now if I forget what the layer number is, but remember that it had the word "Wall" in the name I can just
sealay wall
and like magic I am presented with a list of options which when selected issues a lay= command.

  Likewise if I forget a function key and don't want the entire funkey command list I can just...
seafk odot
and all the appropriate function keys are presented.

It must be nice to have a limitless unhindered storage capacity, but for those of us who have more useful years in the past than in the future, search is our friend.

They both essentially work the same, just one searches the layer names and one searches the function keys.  Here is an example but you get the idea.

print 'seafky.py modified 9:11 AM 7/27/2017'# Display FKeys matching search string
'''
Copyright 2017 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.

Prompts user for Fkey name (or any part) then displays Fkey
'''

def lines2lists(AListOfDataLines):
    DataList=[]
    for Line in AListOfDataLines:
        DataList.append(Line.split())
    return DataList

PrintList=[]

if VrArgs:
    SearchText=VrArgs[0]
else:
    SearchText=PyVrGui().InputDialog ('Search String', 'Partial Funkey Name')[1]

FunKeyFile = open(VrCfg().GetFkeyFileName() , 'r')
FunkeyData = FunKeyFile.readlines()
FunKeyFile.close()
DataList = lines2lists(FunkeyData)

for DataLine in DataList :
    if len(DataLine) == 3 :
        if DataLine[1]=='KeyName':
            if DataLine[2].upper().count(SearchText.upper()):
                PrintList.append(DataLine[2])

if PrintList:
    PromBox = VrPromBox ("Run Fkey", 60, 1)
    PromBox.AddList ("Function Key", 40, len(PrintList)+1, 0)
    for FkeyItem in PrintList:
        PromBox.AddListItem (FkeyItem)
    if (PromBox.Display(0) == 0):
        RunFkey = PromBox.GetListByPrompt ("Function Key")
    else: RunFkey=0
    if RunFkey : PyVrGui().PushKeyin(RunFkey)
    else: PyVrGui().MsgBox('No Match to {:s}'.format(SearchText), 'Matching Layers')

Going from Zero to "How did I ever get along before that" in three minutes

The other day I heard about a little indicator app ( aka bookmarks-indicator ) for Ubuntu that sits nicely on the task bar and allows you to drill down through the disk hierarchy by simply mousing over the icon.  It displays the folder names as you move through the system, auto expanding as you pause over a file name, showing individual files in each folder and opening them when you click.  One minute I didn't even know I wanted this ability and three minutes later I wondered how I ever survived without it.  The fact that is was just a couple dozen easy to read lines of python was even funner.  

Jump ahead to this morning and I came across a bunch of lines that just needed to have one particular point removed from them.  Long story, but on this particular job it happens on a fairly regular basis.  Well I've always had a program to straighten a line, removing all points between to identified locations, but I just needed to delete one point, and doing it with a single click would just be that much easier.  Granted EditVertex set to delete is 95% of the way to what I want but honestly sometimes you just need to step outside for a break to grab that last 5% especially if you can do so with a few lines of code.  And let's be honest it's funner than working anyway.

Hence dellpt.py

print 'dellpt.py modified 9:06 AM 9/20/2016'
# Delete Line Point
'''
Copyright 2016 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.

Deletes a single point identified on a line.
'''

Ws=PyVrWs()
Line=PyVrLine()
while Line.Id() != -1:
    Ws.UndoBegin(WsNum,"dellpt")
    PntNum=Line.GetIdPoint ()
    Line.DelPoint (PntNum)
    Line.ReRec()
    Line.Plot()
    PyVrGr().Replot()
    Ws.UndoEnd(WsNum)

Simple tools for simple tasks (an appreciation for....)

Maybe because I am fairly simple minded, I have always had a deep appreciation for "The Unix Philosophy" which to my thinking involves the idea that if you just need to solve a simple little problem quickly you just throw together a simple little tool to do it but at the same time make sure that there is a way to chain these little tools together in a way that can make the aggregate a truly useful and maybe powerful solution. The integration of this idea into the programs I interact with minute by minute all day every day (and the fact that the first versions were on Unix) make Mike Kitaif one of my heroes.  Yes almost every function, option, toolset has an easy, GUI, comprehensive way of dealing with it but at the same time almost all have a simple command equivalent with arguments that can be built upon using function keys, macros, or incredibly the python interpreter. When you combine this with the fact that rather than a monolithic database for hanging on to all the settings they are usually stored in tiny easily accessible (sometimes even text) files you get the opportunity to make the system do almost anything you can dream.

A great example of all this came up recently when I realized that under certain conditions I wanted to under certain conditions be able to view a particular vertical slice of LiDAR data and set the coloring based on elevation sub-slices.  There is of course a way to do this using the GUI for the poidis command.

There are also fairly straightforward ways to set all these things from the command line.
http://www.vrmapping.net/help5/index.html?point_display.htm

However, what I really wanted was the ability to input a minimal amount of information (base elevation and total slice thickness) and then have it display exactly what I wanted with the colors set to the sub slices (thickness divided by total number of available slices, or 20 currently) and here is where the python comes in. Granted the programming is about at the level of an 8 year old with a Raspberry Pi but you get the idea.

print 'viewsl.py modified 6:06 AM 5/6/2016'
# Views a slice of lidar with color by slice elevation.
'''
Copyright 2016 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.

A little more user friendly option to the Vr alternative which uses an
external parameter file which can be referenced by other programs
like bumvsl

'''
import os
import man
Ws=PyVrWs()
WsNum=Ws.Aws()
Gui=PyVrGui()

MinZ = Ws.GetMinMax (WsNum)[2]
MaxZ = Ws.GetMinMax (WsNum)[5]
AverageZ = (MaxZ+MinZ)/2.0

#In case the parameter file was never created before or can't be read for some reason
BaseZ=AverageZ
DelZMin=MinZ-AverageZ
DelZMax=MaxZ-AverageZ
Thickness=MaxZ-MinZ

ParFileName=VrCfg().GetVrHomeDir()+'/hostdir/viewslice.par'
if not os.path.isfile(ParFileName):
     open(ParFileName,'w').write('{:.2f} {:.2f} {:.2f} {:.2f}'.format(BaseZ,DelZMin,DelZMax,Thickness))
elif not os.path.getsize(ParFileName):
     open(ParFileName,'w').write('{:.2f} {:.2f} {:.2f} {:.2f}'.format(BaseZ,DelZMin,DelZMax,Thickness))
Params=open(ParFileName,'r').read().split()
BaseZ=float(Params[0])
Thickness=float(Params[3])

#In case the user passes in the numbers from the command line.
if VrArgs:
      BaseZ=float(VrArgs[0])
      Thickness=float(VrArgs[1])

PromBox = VrPromBox ("Set Z Slice", 30, 1)
PromBox.SetFocus()
PromBox.AddDouble('BaseZ',BaseZ,2)
PromBox.AddDouble('Thickness',Thickness,2)

if not VrArgs:
      if PromBox.Display(1)==0:
          BaseZ=PromBox.GetDoubleByPrompt('BaseZ')
          Thickness=PromBox.GetDoubleByPrompt('Thickness')
open(ParFileName,'w').write('{:.2f} {:.2f} {:.2f} {:.2f}'.format(BaseZ,DelZMin,DelZMax,Thickness))

#Just use the numbers to pass along arguments to the built in commands.
print BaseZ,DelZMin,DelZMax,Thickness
Gui.PushKeyin ('PoiFilZsl  On')
Gui.PushKeyin ('PoiFilZsl  BasZ {:.2f}'.format(BaseZ))
Gui.PushKeyin ('PoiFilZsl  MinDel {:.2f}'.format(0))
Gui.PushKeyin ('PoiFilZsl  MaxDel {:.2f}'.format(Thickness))

NumberOfSlices=20
SliceThickness = Thickness*1.08/NumberOfSlices #1.08 because it minimizes the top white slice.
print MinZ,MaxZ,SliceThickness,NumberOfSlices
Gui.PushKeyin ('PoiColSli  On')
Gui.PushKeyin ('PoiColSli  ResZ')
Gui.PushKeyin ('PoiColSli SliZ {:d} {:.2f} {:.2f}'.format(1,0,BaseZ))
CurrentElevation=BaseZ
for CurrentSliceNumber in range(2,NumberOfSlices):
    Gui.PushKeyin ('PoiColSli SliZ {:d} {:.2f} {:.2f}'.format(CurrentSliceNumber,CurrentElevation,CurrentElevation+SliceThickness))
    CurrentElevation+=SliceThickness
Gui.PushKeyin ('PoiColSli SliZ {:d} {:.2f} {:.2f}'.format(NumberOfSlices,CurrentElevation,1000000))

PyVrGr().Replot()

*Read the Wikipedia article if for no other reason than to better appreciate some of the men who helped make this all possible.

30 Second Functions

When is it faster to write the function than to even find out how to do it in the help.  Well ok, I did actually search and didn't find out how, but here is what would have been faster.

I use handwheels,  I wanted an easy way to switch the z input from the footdisk to the right handwheel and back.  There is probably a way to pass an argument into a macro but I couldn't find it so enter....

zhand.py

print 'z handwheel'
VrCfg().SetHwLocZ (2)

The beauty of VrCfg is that it modifies the environment in real time just like calling any of the functions that would do the same thing via a dialog box.

Saying anything else would just not have a point, except, always be thinking of ways to extend the usefulness of existing processes.

I Get By With A Little Help From My Friends

Talking to a friend the other day it occurred to me how fantastic it would be if Vr could import and display OpenStreetMap data.  Nothing fancy mind you but maybe create an image that you could use as a background for a quick check of the surrounding area or even flight prep.  Wow that sounds like a lot of work unless of course someone else has already done it.  One of the many beauties of Python is the number of fantastic libraries out there if you just spend a little time tracking them down.  If you take that a step farther and are willing to pass a little of the responsibility outside python using system calls or batch files the possibilities are nearly endless.  Let's just say I can now display an OpenStreetMap map image in Vr and it only took about 90 lines of python.  Here are a few of the major helpers besides the map site which is so worthy of support, contribution, and praise.

I know I have mentioned it before but I can't say enough about pyproj which I use to convert my local coordinates to geographic for sending out to various places.

Once I have the corners in Lat / Long the rest will be done by a batch file which is about the only thing my python script creates.  It could be done in other ways but I liked the idea of creating something that could stick around and be easily modified for use again.  So the batch file next makes a call to the openstreetmap.org API for grabbing the map data.

Once the OSM file is in house it will be converted to a bitmap using Maperitive.  I only use such a tiny fraction of the capabilities of this amazing program, but for my purposes using it's command line version with a custom generated script will fit the bill perfectly.

The original bitmap is created in Web Mercator and can also throw in a TFW and KML file for generic use or display by anybody if you want to send them along.  I however want to display it behind my map in a local coordinate system in whatever system I choose so the last step is to re-project the image using GDAL.

Is this the most advanced, professional, elegant way of doing it?  Of course not, you obviously don't know me very well.  As usual it is something you can hack together in just a short time with minimal duplication of effort that gets the job done.  Hooray for real programmers!!!

Don't forget to support your favorite projects with contributions or cash when possible.

Note to self: Step backward through LiDAR when removing, and repr ( round ( a number )) is awsome

First, apparently I never noticed this before but while trying to step though a bunch of LiDAR points to remove a certain class it apparently changes something dynamically (like the point number in the buffer).  I noticed this because I was doing a

    for PointNum in range(Punt.GetCount()):
        if Punt.Cla(PointNum) == Whatever:
            Punt.DelPunt(PointNum)

To get rid of a certain type of point, in this case based on class.  What I noticed though was that each time it ran there were a bunch of points left over (half??).  I decided that even though I wasn't doing a record until the loop was done, the point number in the buffer must change as points are deleted.  So what to do?  How about starting at the last point and removing them from the top end of the buffer so as you work backwards the previous point numbers would not have changed.  I came up with

for PointNum in range(Punt.GetCount(),-1,-1):

and it seems to work fine.

Second and totally unrelated I needed to rebuild a tile naming structure for a local dataset.  There is a tile scheme that is based on 5000 foot tiles that is easy enough to get the names of.  The problem is that there is also a subset of 1250 foot tiles that starts from the corner of each 5k block but the names are based on 1000 foot increments.  When you digitize a random point in a block you can't necessarily round down to the lower 1000 because they just truncate the thousands from a rounded 1250.  Ok hard to imagine but the corners would be 0, 1250, 2500, 3750 so a number 3749 would fall in a 2500 block and truncate to 2000 meaning the tile would have 2 in the name, not to mention dropping a million off one of the coordinates.  Now for the why I love python part!

repr(round(1773691/1250)*1250)[1:4]

yields '772' just like is necessary to build the proper tile name. 

Note to self: Find a reason to use geopy

Full disclosure; This has nothing to do with anything except that I was playing with a really cool module today and need to make sure I remember what I just learned.  The module is geopy (which by the way I installed using pip for the first time ever ) which can be used to geocode natural language addresses using a variety of APIs.  Don't know when, how, or why I'll use it but I'm certain I will.

Info can be found at https://pypi.python.org/pypi/geopy
Doc at http://geopy.readthedocs.org/en/latest/

Here is what I just did (using my current address)

>>> import geopy
>>> geopy.geocoders.GoogleV3().geocode("4090 weaver ct s 43026")
Location((40.035917, -83.144145, 0.0))
>>> geopy.geocoders.GoogleV3().geocode("4090 weaver ct s 43026").raw
{u'geometry': {u'location': {u'lat': 40.035917, u'lng': -83.144145}, u'viewport': {u'northeast': {u'lat': 40.03726598029149, u'lng': -83.1427960197085}, u'southwest': {u'lat': 40.03456801970849, u'lng': -83.1454939802915}}, u'location_type': u'ROOFTOP'}, u'formatted_address': u'4090 Weaver Court, Hilliard, OH 43026, USA', u'place_id': u'ChIJpUOLF-iTOIgRT-Vqu5rHhLo', u'address_components': [{u'long_name': u'4090', u'types': [u'street_number'], u'short_name': u'4090'}, {u'long_name': u'Weaver Court', u'types': [u'route'], u'short_name': u'Weaver Ct'}, {u'long_name': u'Northwest Industrial Complex', u'types': [u'neighborhood', u'political'], u'short_name': u'Northwest Industrial Complex'}, {u'long_name': u'Hilliard', u'types': [u'locality', u'political'], u'short_name': u'Hilliard'}, {u'long_name': u'Norwich', u'types': [u'administrative_area_level_3', u'political'], u'short_name': u'Norwich'}, {u'long_name': u'Franklin County', u'types': [u'administrative_area_level_2', u'political'], u'short_name': u'Franklin County'}, {u'long_name': u'Ohio', u'types': [u'administrative_area_level_1', u'political'], u'short_name': u'OH'}, {u'long_name': u'United States', u'types': [u'country', u'political'], u'short_name': u'US'}, {u'long_name': u'43026', u'types': [u'postal_code'], u'short_name': u'43026'}, {u'long_name': u'1197', u'types': [u'postal_code_suffix'], u'short_name': u'1197'}], u'partial_match': True, u'types': [u'street_address']}

LMSTFY ( Let Me Sealay.py That For You)

I should have mentioned this long ago before you invested time in reading my blog, but here goes... Fact is I'm not particularly bright, I'm just inquisitive (and as mentioned previously constructively lazy).  This means a couple of things; first I love to find an easier way to do almost anything, and second sometimes I spend more time researching than I would if I had just done the task (oh but next time it will be so much faster). Because of this, sometimes I know things, and when you know things people like to ask you things because you might know that particular thing they want.  Yes almost daily I am tempted to just say "Let me Google that for you" because that would be much easier for me than actually explaining it (reference lazy thing above).

I have about 150 out of the 240 layer numbers that I use daily in mapping projects memorized.  Over time though you need to add new layers or break them down to a more granular level so consequently you can see that there are a bunch I don't know.  I could memorize them (and a few I will) but the majority of the others I just want an easy way to find the rare one I need.  Hence sealay.py or "Search Layer Names".

print 'sealay.py modified 7:31 AM 4/16/2015'
# Display layers matching search string
'''
Copyright 2015 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.

Prompts user for layer name (or any part) then displays layer
numbers and full names of layers matching string.
'''

def lines2lists(AListOfDataLines):
    '''
    Function readlines returns an entire file with each line as a string in a
    list of data.  This function will convert each string into a list of words,
    then return a list of lists. Why? Just because I like to work this way.
    Example:            lines2lists(['first line','the second line','or sentences'])
    Would return:       [['first','line'],['the','second','line'],['or','sentences]]
    '''
    DataList=[]
    for Line in AListOfDataLines:
        DataList.append(Line.split())
    return DataList

PrintString=''
PrintList=[]

#If sent argument use as search string else prompt.
if VrArgs:
    SearchText=VrArgs[0]
else:
    SearchText=PyVrGui().InputDialog ('Search String', 'Layer numbers by name')[1]
   
LayerFileName = VrCfg().GetLayerNameFile ()
LayerFile = open(LayerFileName , 'r')
LayerData = LayerFile.readlines()
LayerFile.close()
LayerList = man.lines2lists(LayerData)

#Search through data, if string found add it to a multi-line printable string.
for DataLine in LayerList :
    if len(DataLine) > 1 :
        if DataLine[1].count(SearchText.upper()):
            PrintString=PrintString+DataLine[0]+' '+DataLine[1]+'\n'
            PrintList.append(DataLine[0]+' '+DataLine[1])

# Two possibilities, the first commented line just displays the found layers
# The rest will let the user select one correct entry and send a LAY= command.
if PrintString:
#    PyVrGui().MsgBox(PrintString , 'Matching Layers')
    PromBox = VrPromBox ("Set Layer", 30, 1)
    PromBox.AddList ("Layer", 40, len(PrintList)+1, 0)
    for LayerItem in PrintList:
        PromBox.AddListItem (LayerItem)
    if (PromBox.Display(0) == 0):
        SetToLayer= PromBox.GetListByPrompt ("Layer")
    if SetToLayer : PyVrGui().PushKeyin('lay={:s}'.format( SetToLayer.split()[0]))
else :
    PyVrGui().MsgBox('No Match to {:s}'.format(SearchText), 'Matching Layers')

You gotta love a good library

The Columbus Metropolitan Library is arguably one of the best libraries in the nation.  I use it so much I have had my library card memorized for over 25 years.  The idea that there is information and services there at my fingertips any time I want is fantastic.  I was thinking about this yesterday when I started pulling functions out of PyVrGeom().  Libraries (modules) can be everything from little bits of code snippets that you want to reuse over and over to complex collections of classes that you could only dream up but never realistically code yourself. 

For example PyVrGeom is just a collection of math that some would find easy enough to code, but why do so when somebody has done the heavy lifting and offers you the calls to do it with some preexisting code.  In this case very professionally created and presented, but it is just as easy to hack together some snippets that just save you typing.  For example when I was working on a KML export routine I didn't need to learn something complex and universal (which libs probably exist), I just knew what it looked like and wanted to replicate those entities by passing along a few basic parameters.  In this case I slapped together some code to create the KML entities, but at the same time the coordinates needed to be in geodetic coordinates so I had to grab a projection library that someone way sharper than me had already put together. Then it struck me that a KMZ is just a file zipped using the standard algorithms so if I import that one I can just as easily write both KML and KMZ files.

I'm not going to offer any code this time just another reminder that there are some good reasons to install python if you are going to do some scripting.  Here are some of the libraries that I really enjoy and use often.

pyproj - for converting local coordinates to global (or from one to another).
liblas - for all things LiDAR related.
zipfile - for reading and writing compressed files.
xlrd/xlwt - for directly reading and writing Excel spreadsheets.

Along with some of the standard system libraries that you can find a million uses for like sys, os, time, and math.

And here is a tip I'll add for free, check out http://www.lfd.uci.edu/~gohlke/pythonlibs/

It's all part of the process (or subprocess as the case may be)

On occasion I make the case for installing python in order to extend the usefulness of the environment in Vr, or at least getting hold of some of the most useful standard libraries.  Well here is another good reason.  This is just a snippet that I will use later but I figure this is as good a place as any to keep track of it.

Before too long I want to write a script that will depend on running an external process, then interacting with the results of the process (no spoilers yet).  In order to do that I'm going to make something similar to the following call to process some files leaving the results in the working directory. The research I did here was to make sure I had a way to not proceed until the original process was complete.  Later I'll need to make sure I can run a shell command and capture the resulting output but this will do for now.

import subprocess
process = subprocess.Popen('lasboundary -i *rgb1.las -otxt')
process.wait()
print process.returncode
print '\a'

the last line is just to beep so I knew when it was done.

Reminder snippets - liblas and reading laz files

The other thing that I used to do with this blog is post things that I wanted to remember.  At one point the Vr forums were a good place to do this because it also gave an opportunity to share and collaborate.  Since that doesn't really seem to be happening any more I suppose I could drop things here and then come back to them later if necessary.  Yes I could just do this in a document but who knows maybe a new opportunity for some kind of sharing community will pop up again.

In this case I just want to remind myself that if Vr never adopts the .LAZ LiDAR compression format for reading, it would still be possible using liblas. This is simple a little interactive session that shows it is theoretically possible using my own data.

There may indeed be other better ways, but this seems to work fine.

>>> import liblas
>>> f=liblas.file.File(r'c:\tmp\N1915250.laz',mode='r')
>>> f

>>> p=f.read(100)
>>> p.classification
5
>>> p.x
1919853.35
>>> f.header.compressed
True
>>> f.filename
'c:\\tmp\\N1915250.laz'

5 Minute Functions - 2

I haven't stopped coding, I've just stopped talking about it.  Again more a function of the fact that I doubt anybody really cares.  I noticed the other day that I have slightly over 300 python programs that have accumulated over the years.  Certainly there are the ones I am really proud of like the ones that translate a Vr file into a KML file or a TIN into  LandXML format because of the crucially useful functionality they add.  Then there are the ones I use multiple times an hour, like double points on a line, drive to points by user specified parameters, storing helpful window states, adding clamping points to DTM lines, and on it goes.  Lately I have noticed that I also just stop sometimes and write a quick program because as I mentioned in a previous post of a similar title, I can spend 5 minutes programming and save 15 minutes (and lots of wear and tear on the mouse clicking finger).  Not only is it worth it but it is also fun, keeps my brain lubricated, and gives me a base in case I find the task useful and want to come back later to add options, dialogs, or arguments that would make it more universal.

Here is an example, I have a few thousand points that got dumped into a file without regard to any distinction.  about 300 of them would be really helpful if they were in a certain layer and I don't want to miss any of them.  The only thing that makes them different is that there is a common string of text inside the feature code.  There may very well be an easy way to do this with Vr, but after a couple of tries I couldn't come up with a parameter that worked so I thought, "why not just write a quick script?".  In this instance I'll just hard code the search string and resulting changes. If I find it useful I may come back later and add a simple dialog to make it more agile.

print 'fcglob.py modified 6:23 AM 3/11/2015'
# Globally change something based on text in it's feature code
'''
Copyright 2015 Dennis Shimer
No warranties as to safety or usability expressed or implied.
Free to use, copy, modify, distribute with author credit.

Simple beginning of a function to change entities based on their
feature code.

Variables of interest:
    None
'''
Ws=PyVrWs()
Sym=PyVrSym()
WsNum=Ws.Aws()
Ws.UndoBegin(WsNum,'fcglob')
for EntNum in range ( Ws.GetSymCount(WsNum)):
    Sym.Load (WsNum, EntNum)
    if Sym.GetFc().lower().count(' fh'):
        print EntNum
        Sym.SetLayer(1042)
        Sym.SetGpoint(21)
        Sym.ReRec()
print ('\a')
Ws.UndoEnd(WsNum)

 It started out as about a half dozen lines in pyedi and doesn't amount to much more now. Nothing fancy, just something I can come back to if I need it again, and in this case it really helped.

I don't know, maybe I'll just start dumping more of these here. Doesn't hurt anything and though most of them aren't universally interesting it might give somebody an idea.

Lazy - Efficient, Potato - Potaaahto

I don't really tend to share my most complex and production related scripts because this is more of a tutorial and tip related blog.  In that same vane I have mentioned in the past that there are times that you do the same task several times or come up with a new procedure that you implement in several jobs that just begs for a simple script.  This is one of those times.

I have found over the years that having a simple command that saving the current states of things so you can recall them later is a big time saver.  For example save the current window display parameters so you can come back to the same settings later, or save the layer on/off states for easy recall.  There are ways to do stuff like this but I just want one two letter macro to save something and another to restore it.  Lately I have been working with 5-6 files at a time and if I jump out of the program or switch from VrOne to VrTwo I want to open the same files.  Hence StoVr.py, and GetVr.py to take a snapshot of open files and re-open them if I need to.

print 'stovr.py modified 9:08 AM Friday, August 30, 2013'
# Store Vr filenames of open files for later recall
'''
Copyright 2013 Dennis Shimer
No warranties as to safety or usability expressed or implied.
Free to use, copy, modify, distribute with author credit.

Echos the names  of currently open files to the command windows
and saves them for later recall.
'''

Ws=PyVrWs()
WorkSpaceCount=Ws.GetWsCount()
print 'File names saved'
if WorkSpaceCount:
    SaveFileNamesFile=open(VrCfg().GetVrHomeDir()+'\\hostdir\\SavedVrFileNames.txt','w')
    for WsNum in range(WorkSpaceCount):
        SaveFileNamesFile.write(Ws.GetFileName (WsNum)+'\n')
        print Ws.GetFileName (WsNum)
    SaveFileNamesFile.close()
else: print 'No workspaces open'

print 'getvr.py modified 9:08 AM Friday, August 30, 2013'
# Get Vr filenames previously stored
'''
Copyright 2013 Dennis Shimer
No warranties as to safety or usability expressed or implied.
Free to use, copy, modify, distribute with author credit.

Echos the previously saved names to the command window with
workspace numbers as they are opened.
'''

Ws=PyVrWs()
SaveFileNamesFile=open(VrCfg().GetVrHomeDir()+'\\hostdir\\SavedVrFileNames.txt','r')
FileNameDataList=SaveFileNamesFile.readlines()
for FileName in (FileNameDataList):
    PyVrGui().PushKeyin ('opevr '+FileName.rstrip())
    print Ws.GetWsCount(),FileName.rstrip()
SaveFileNamesFile.close()

Let me look that up in the dictionary for you.

There have been some changes to the VrPunt methods over time.  I have suspected that not everything is implemented exactly the way it is documented, and sometimes a name could be different or missing.  I have discovered though that the real magic is just in working with the dictionary which stores all the relevant LiDAR point attributes.  The LiDAR data provider that we most often use stores the individual flight line tag in the "SOU" tag, Vr however will only filter by the "FLT" tag.  Since FLT isn't used in our data, changing it to reflect the SOU value doesn't seem to have a downside.  Below is a simple program that runs through all the point data and copies SOU to FLT.  Now I can use the flight line filter to turn the points on and off.

print 'souflt.py modified 12:44 PM 8/27/2013'
# Copy LiDAR Source attribute to Flight attribute
'''
Copyright 2013 Dennis Shimer
No warranties as to safety or usability expressed or implied.
Free to use, copy, modify, distribute with author credit.

Reads the "Source" in a PuntA dictionary and sets the "Flight" to the same value.

'''

Ws=PyVrWs()
Punt=PyVrPunt()
Gui=PyVrGui()
WsNum=Ws.Aws()
PointBufferCount=Ws.GetPuntBufCount(WsNum)

Ws.UndoBegin(WsNum,"Source2Flight")
Gui.ProgInit('Points',Ws.GetPuntBufCount(WsNum))
for PointBufferNumber in range(PointBufferCount):
    Punt.Load(WsNum,PointBufferNumber)
    Gui.ProgSet(PointBufferNumber)
    for PointNum in range(Punt.GetCount()-1,-1,-1):
        PuntA=Punt.GetPuntA(PointNum)
        PuntA['Flt']=Punt.Sou(PointNum)
        Punt.ChgPuntA (PointNum, PuntA)
    Punt.ReRec()
Gui.ProgReset()
PyVrGr().Replot()
Ws.UndoEnd(WsNum)

 So the key is to load up the PuntA to make sure it is populated with the correct values, then modify the particular attribute and save it back to PuntA. When all the points are done re-record the entire Punt buffer.  As of this writting the dictionary is organized as follows.

"Lay"     = 1   - Layer                    (1-30001)
"Int"     = 0   - Intensity                (0-65535)
"Dsp"     = 1   - Display flag             (0-1)
"Ret"     = 0   - Return number            (1-5)
"Nre"     = 0   - Number of returns        (1-5)
"Sdf"     = 0   - Scan direction flag      (0-1)
"Edg"     = 0   - Edge of flight line flag (0-1)
"Cla"     = 0   - Classification           (0-31)
"Syn"     = 0   - Synthetic flag           (0-1)
"Key"     = 0   - Key-point flag           (0-1)
"Del"     = 0   - Delete flag              (0-1)
"Ang"     = 0   - Scan angle               (-90 - +90)
"Flt"     = 0   - Flight number            (0-255)
"Sou"     = 0   - Point source Id          (0-65535)
"Red"     = 255 - Red component color      (0-255)
"Green"   = 255 - Green component color    (0-255)
"Blue"    = 255 - Blue component color     (0-255)
"GpsTime" = 0.0 - GPS Time                 (Double precision number)

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 showll.py 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

preserve_units=True 

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
 Gui.DspMsg0(repr(Projection(x,y,inverse=True)))
 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)

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

try: 
 MyApp.SetDigCB (MyAppDigCB)
 Mk.SetKeyCB (MyAppKeyCB)
 MyApp.KeepRunning ()
 
 Mk.Delete ()
except:
 print "Exception during application:", sys.exc_info()[0]
 Mk.Delete()
 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)

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/

We all have our limits

Note that I am trying a new method of posting code, hope it works better but of course the original code itself can always be downloaded from the website where it is stored.

Quite often when writing some sort of import function I would like to limit the data imported to certain physical constraints. I started out asking the user if they want to only act on coordinate values that fall within a bounding box, then if they said yes, prompt them to ID the bounding line. This works great as far as it goes, as part of the ID function I would set a flag (in this case BoundrySet) to 1 then later in the program if BoundrySet was true the IsPointInside method of the Line class could be called. Something like....

if BoundrySet :
if BoundingLine.IsPointInside (X,Y): Do something

But then I got thinking, there are several easy ways to check for bounding areas without having a line already digitized. How about using the screen boundaries, or digitizing a rectangle based on the lower left and upper right corners, or for that matter just start digitizing points and create a line in memory on the fly. Well it turns out that I already had the snippets I needed in other functions so just copied them out, cleaned them up a bit (and I do mean a bit, it still isn't particularly pretty, but it works) and created a basic function that I can paste into any program where spatial filtering is helpful.

There are a couple of caveats. Using the rectangle select mode requires that you select the lower left point first, then the upper right based on the screen orientation. I'll probably clean this up some day and allow for more flexibility but it falls into the "quick hack for some purpose, and I know how it acts" category of code that gets the job done but isn't as attractive as it should be. Using the "Digitize" mode allows the user to start clicking on positions and builds the boundary line on the fly, if it would be helpful it would be just as easy to record the line for future interaction but for now I just throw it on the screen. The nice thing about this dialog is that it behaves just the way you expect so if the focus is in the drop down and you hit an 'r' rectangle mode will be activated and so on.

print 'SetBoundary.py modified 8:00 AM 5/15/2009'
# Set spatial bounding for processing
'''
Copyright 2009 Dennis Shimer, M.A.N. Mapping Services Inc.
No warranties as to safety or usability expressed or implied.
Free to use, copy, modify, distribute with author credit.

Really just a container for the code snippet that would do spatial
filtering as part of a larger program. Just written as working code
here for testing purposes.

Variables of interest:
None
'''
import math

BoundingLine=PyVrLine()
Gr=PyVrGr ()
Gui=PyVrGui()
BoundrySet=0

PromBox = VrPromBox ("Set Boundry", 30, 1)
PromBox.AddCombo ('Boundry', 5, 0, 0)
PromBox.AddComboItem ('Line')
PromBox.AddComboItem ('Screen')
PromBox.AddComboItem ('Rectangle')
PromBox.AddComboItem ('Digitize Points')
PromBox.AddComboItem ('None (whole file)')
PromBox.SetFocus()
if PromBox.Display(1)==0:
BoundryType=PromBox.GetComboByPrompt ('Boundry')
if BoundryType==0:
   BoundingLine.Id()
   BoundrySet=1
elif BoundryType==1:
   for corner in range(4):
       x,y,z=Gr.GetWinCorner (corner,0)
       BoundingLine.AddPoint(x,y,z)
   BoundrySet=1
elif BoundryType==2:
   Stat,llx,lly,z=Gui.GetCoord('Boundry')
   Gr.DrawMode (1,-1)
   Gr.PenNum (1)
   Gr.DrawMarker (llx,lly,z,MARK_CROSS,Gr.GetWinScale ()*.2,-1)
   Stat,urx,ury,z=Gui.GetCoord('Boundry')
   Gr.EraseMarker ()
   ScreenAngle=Gr.GetWinRot()[2]
   Angle=math.atan((ury-lly)/(urx-llx))-Gr.GetWinRot()[2]
   Diagonal=PyVrGeom().Dist(llx,lly,urx,ury)
   Height=math.sin(Angle)*Diagonal
   if Height < 0.0: Height=Height*-1.0
   Width=math.cos(Angle)*Diagonal
   if Width < 0.0: Width=Width*-1.0
   BoundingLine.AddPoint(llx,lly,z)
   lrx=llx+(Width*math.cos(ScreenAngle))
   lry=lly+(Width*math.sin(ScreenAngle))
   BoundingLine.AddPoint(lrx,lry,z)
   urx=lrx-(Height*math.sin(ScreenAngle))
   ury=lry+(Height*math.cos(ScreenAngle))
   BoundingLine.AddPoint(urx,ury,z)
   ulx=urx-(Width*math.cos(ScreenAngle))
   uly=ury-(Width*math.sin(ScreenAngle))
   BoundingLine.AddPoint(ulx,uly,z)
   BoundingLine.AddPoint(llx,lly,z)
   BoundingLine.Plot()
   BoundrySet=1
elif BoundryType==3:
   Stat=0
   while Stat==0:
       Stat,x,y,z=PyVrGui().GetCoord('Dig Boundry')
       BoundingLine.AddPoint (x,y,0.0)
       BoundingLine.Plot()
   BoundingLine.Close(2)
   BoundingLine.Plot()
   BoundrySet=1