Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

=== SUMMARY === I have created a process to import a LAS point cloud, clip it, select the ground points and save the output as an XYZ.txt for clients to use. This uses a python based script (will become a tool) that utilizes some of the better features in SAGA GIS as well as arcpy.

SAGA seems to first create/save the spc (saga point cloud) and then re-import it before creating a shape file. Can I get it to work only in the memory?

Also once the spc is cut by the AOI, I want it to split the file into 200,000 record segments (max allowable by some clients sofware) -can this be done?

=== UPDATE 2 === based on @johanvdw posting of Wiki with answers ==

I'm still having issues with what the sys arg values are and how to set the input/output. I now have the following test code...I have added questions to it.

import saga_api, sys, os, arcpy

##########################################
def las2spc(fLas):
# fLas: python string containing the filename
# returns a pointer to a pointcloud
    print saga_api.SAGA_API_Get_Version()

    Library = saga_api.CSG_Module_Library()
    if not Library.Create(saga_api.CSG_String('C:\SAGA\modules\io_shapes_las.dll')):
        print 'unable to load SAGA module library'
        return 0

    Module = Library.Get_Module('Import LAS Files')
    print "Module Loaded"
    pSpc  = saga_api.SG_Create_PointCloud()
    Parms = Module.Get_Parameters() # default parameter list
    Parms('FILE').Set_Value(saga_api.CSG_String(fLas)) #filename in saga string format
    Parms('POINTS').Set_Value(pSpc) #pointcloud data object as required by the module (see above)
    Parms('T').Set_Value(0) #This are booleans, either on or off
    Parms('i').Set_Value(0) #T eg is the time, see above
    Parms('a').Set_Value(0)
    Parms('r').Set_Value(0)
    Parms('c').Set_Value(0)
    Parms('u').Set_Value(0)
    Parms('n').Set_Value(0)
    Parms('R').Set_Value(0)
    Parms('G').Set_Value(0)
    Parms('B').Set_Value(0)
    Parms('e').Set_Value(0)
    Parms('d').Set_Value(0)
    Parms('p').Set_Value(0)
    Parms('C').Set_Value(0)
    Parms('VALID').Set_Value(0)
    Parms('RGB_RANGE').Set_Value(0)

    if not Module.Execute():
        print 'module execution failed'
        return 0
    print 'module successfully executed'

    return pSpc 


if __name__ == '__main__':

# I have commented the following out as it wasn't clear what it does. Are the sys.arg values from running this as a tool?

##    if len( sys.argv ) != 2:
##        print 'Usage: las2shape.py <in: shape> <out: pointcloud>'
##    else:
        fLas   = '335000_8084000_1K_class.las' 

I put in the name of a las file in the same directory as the py file.**

        if os.path.split(fLas)[0] == '':
            fLas   = './' + fLas
        pSpc = las2spc(fLas)
        fSpc   = 'file.spc' 

this is the name I want to give the output

        if os.path.split(fLas)[0] == '':
            fSpc   = './' + fSpc
            print "Saving"
        pSpc.Save(saga_api.CSG_String(fSpc))
        print "Saved"
#or instead of saving: use pSpc as on of the inputs to your next module; the wiki gives an example.

which gives the following result.

SAGA Application Programming Interface - Version: 2.0.7 Module Loaded module execution failed Traceback (most recent call last): File "P:\2011\Job_193_LAS_2_SHP_Procedure\Working\NewScript\SAGA.py", line 57, in pSpc.Save(saga_api.CSG_String(fSpc)) AttributeError: 'int' object has no attribute 'Save'

=== UPDATE 1 === based on @johanvdw posting of Wiki with answers ==

Based on the input that we can't save in_memory using the cmd line, I have switched to python. I am still having issues with figuring out what some of the code means.

  1. What do all the sys.arg[1-4] values mean? Is this designed to run as an arcgis tool?
  2. "OUTPUT .Save(saga_api.CSG_String(sys.argv2))" is the last line -

When I try the code below it gives me an error -

# Python script template for SAGA module execution (automatically created, experimental)

import saga_api, sys, os, csv

def SAGA_importLas(SHAPES, ZFIELD, OUTPUT):
    print saga_api.SAGA_API_Get_Version()

    Library = saga_api.CSG_Module_Library()
    if Library.Create(saga_api.CSG_String('C:\SAGA_win32\modules\pointcloud_tools.dll')) == 0:
        print 'unable to load SAGA module library'
        return 0

    Module = Library.Get_Module('Point Cloud from Shapes')
    #print "Opened module"
    Parms = Module.Get_Parameters() # default parameter list
    Parms('SHAPES').Set_Value(0) # input NOT optional shapes
    Parms('ZFIELD').Set_Value(2)
    Parms('OUTPUT').Set_Value(1)

    print SHAPES, ZFIELD, OUTPUT
    if Module.Execute() == 0:
        print 'module execution failed'
        print sys.stdout
        return 0
    print 'module successfully executed'
    return 1

import os, arcpy, sys, datetime

top = os.getcwd() # change to a specific path if required.
# This otherwise starts with the directory the script is in.
RootOutput = top # change if you want output somewhere else
FileTypes=['las'] # add filetypes or * as required.
SearchStrings=['Temp_Pol', 'Spatial_Ex'] # add strings as required.
RootDirectory = r'L:\Vector_Data\LAS'

filecount=0
#successcount=0
#errorcount=0

print "Working in: "+os.getcwd()

List =[]
f = open(RootOutput+'\\Success_LOG.txt', 'a')
f.write("Log of files Succesfully processed. RESULT of process run @:"+\
        str(datetime.datetime.now())+"\n")
f.close()

FileList = csv.reader(open('Request.csv'))
SearchStrings=[]
for File in FileList:
    #SearchStrings.append(File[0]+","+File[1])
    SearchStrings.append(File[0])

for root, dirs, files in os.walk(RootDirectory, topdown=False):
    for fl in files:
      currentFile=os.path.join(root, fl)
      for FileType in FileTypes:
          status= str.endswith(currentFile,FileType)
          if str(status) == 'True':
              for SearchString in SearchStrings:
                  if str(SearchString in currentFile) == 'True':
                    #print str(currentFile)+str(status)       
                    List.append(currentFile)
      filecount=filecount+1

#print 'File List: ' + str(List) + 'containing ' + str(len(List)) + ' files.'

#for item in List:
Item=List[45]
print "Working on : " + Item
if __name__ == '__main__':
##    if len( sys.argv ) != 4:
##        print 'Usage: this_script.py <in: gridfile> <out: gridfile> <in: shapefile> <out: shapefile>'
##
##    else:
        SHAPES = Item
        ZFIELD = ''
        PositionDot = Item.find('_1K_class')
        PositionDir = Item.rfind('\\')
        OUTPUT = Item[(PositionDir+1):PositionDot]+'.spc'


        if SAGA_importLas(SHAPES, ZFIELD, OUTPUT) != 0:
            OUTPUT  .Save(OUTPUT)

Gives result

Working in: P:\2011\Job_193_LAS_2_SHP_Procedure\Working\scripts Working on : L:\Vector_Data\LAS\Rangeview_Tolga\338000_8095000_1K_class.las SAGA Application Programming Interface - Version: 2.0.7 L:\Vector_Data\LAS\Rangeview_Tolga\338000_8095000_1K_class.las 338000_8095000.spc module execution failed

What's happening?

I also don't understand what if "len( sys.argv ) != 4:" - does? I know it's checking if 4 system arguments have been called but why?

Is this in/out grid code that's automatically created wrong? the python of the two process mentioned below also has the same info repeated.

Once it completed it should run Point Cloud Reclassifier / Subset Extractor where the input is the in_memory from the above process which is then used to run Point Cloud from Shape.

=== ORIGINAL POST === new process is only in python.

The py, and bat files are in a zip in google docs.

Thanks for your advise and help,

Code snippets from bat file (I know there is the python saga_api but I can't seem to download it from http://www.saga-project.org/download/python as it gives a 404 error.

@ECHO OFF

SET SAGA_MLB = C:\SAGA\Modules
SET PATH = %PATH%;C:\SAGA

ECHO %1 %2 %3

PAUSE  

C:\SAGA\saga_cmd io_shapes_las "Import LAS Files" -POINTS= -FILE=%3\%1_1K_class.las -T -i -a -r -c -u -n -R -G -B -e -d -p -C -VALID -RGB_RANGE=0

PAUSE
c:\SAGA\saga_cmd pointcloud_tools "Point Cloud to Shapes" -POINTS= -SHAPES=%2\%1_fc.shp 

PAUSE

This is called from a python script that also uses arcpy

[CODE}
import os, datetime
from subprocess import Popen

LASList = ['329000_8119000','330000_8118000','330000_8119000','331000_8118000','331000_8119000']
#LASList = ['329000_8119000']
workingDIR = r'P:\2011\Job_183_Contours_DaneCraze\Working'
outDIR = workingDIR+"\\test3"
AOI=workingDIR+'\\Data_Request_Extents_region.shp'
SourceDIR=r'L:\Data_Admin\temp_dataset_storage\LiDAR\FinalSupply\Atherton_Tolga_Walkamin_Mareeba_Bilwon\LAS'
#SourceDIR=outDIR
spatial_reference = r'C:\Program Files\ArcGIS\Desktop10.0\Coordinate Systems\Projected Coordinate Systems\National Grids\Australia\GDA 1994 MGA Zone 55.prj'



print "Starting process on "+ str(len(LASList)) +" tiles @ "+str(datetime.datetime.now())+ "\n"

for LASTile in LASList:

# SAGA Processes

# Import LAS

    p = Popen("importLAS_exportSHP.bat "+LASTile+" "+outDIR+" "+SourceDIR, cwd=r"P:\2011\Job_183_Contours_DaneCraze\Working")
    stdout, stderr = p.communicate()  
    print "Imported LAS @ "+str(datetime.datetime.now())+ "\n"


### Convert LAS to shp
##
##    p = Popen("PointCloud_shape.bat "+LASTile+" "+outDIR, cwd=r"P:\2011\Job_183_Contours_DaneCraze\Working")
##    stdout, stderr = p.communicate()
##    print "Converted LAS to SHP @ "+str(datetime.datetime.now())+ "\n"
##    

# ARCGIS Processes

### Set local variables
##    in_features = outDIR+"\\"+LASTile+"_fc.shp"
##    clip_features = AOI
##    out_feature_class = outDIR+"\\"+LASTile+"_clip.shp"
##    xy_tolerance = ""
##
### Execute Clip
##    arcpy.Clip_analysis(in_features, clip_features, out_feature_class, xy_tolerance)
##    arcpy.DefineProjection_management(out_feature_class, spatial_reference)
##    print "Executed Clip @ "+str(datetime.datetime.now())+ "\n"
##
### Execute Subset
##    
##    arcpy.MakeFeatureLayer_management(out_feature_class, "lyr")
##    arcpy.SelectLayerByAttribute_management("lyr", "NEW_SELECTION", ' "CLASSIFICA" = 2 ')
##    print "Selected Subset @ "+str(datetime.datetime.now())+ "\n"
##    
##    arcpy.CopyFeatures_management("lyr", LASTile+"_class2")
##    arcpy.DefineProjection_management("lyr", spatial_reference)
##    print "Executed Subset @ "+str(datetime.datetime.now())+ "\n"
##
### Create XYZ - in SAGA
##
##    p = Popen("export_xyz.bat "+LASTile+" "+outDIR, cwd=r"P:\2011\Job_183_Contours_\Working")
##    stdout, stderr = p.communicate()
##    print "Created XYZ txt and exiting @ "+str(datetime.datetime.now())+ "\n"
##

[/CODE}

One reason for me to post this is that I can't convert spc to shp on large clouds as I get the attached error. It runs fine on smaller tiles. enter image description here

share|improve this question
    
I don't have time to give a full answer, but read the wiki, especially the last sentence: you have to edit the file to work. Perhaps it is better to base yourself on the examples in the wiki rather than on the autogenerated file. –  johanvdw Sep 20 '11 at 7:49
    
thanks. I just wasn't sure how much the autogenerated python needed editing. Anyway I re-wrote it and am getting an error on the SAGA execution. Is there a way to see what exactly the error is? Code and result updated above. –  GeorgeC Sep 20 '11 at 23:38
    
I've updated my answer. You are mixing up a lot of things. First of all you are using the wrong module. You are also not setting the parameters to anything sensible. Check the parameters for the module you need as explained in the wiki. –  johanvdw Sep 21 '11 at 7:38
    
I'm still having issues. Mostly as I don't understand - 1) What do all the sys.arg[1-4] values mean? Is this designed to run as an arcgis tool? 2)I/O where do we feed it the input LAS file and the output etc? –  GeorgeC Sep 22 '11 at 3:13
    
google sys.arg, first hit: docs.python.org/library/sys.html ; you run it from the command line, but of course you could also specify the filenames any other way. –  johanvdw Sep 22 '11 at 7:18
show 5 more comments

2 Answers 2

up vote 3 down vote accepted

(I have reworked my answer and included different comments.)

It is not possible to save things to memory when using batch files. Batch files run programs independent of eachother. There is no mechanism to share memory between these programs. It is however possible by using the python interface: in that case you can leave the spc object in memory prior to writing it as a shapefile.

Instead of giving a full answer here, I've updated the SAGA wiki and a short tutorial on creating python scripts, demonstrating how this can be done:

https://sourceforge.net/apps/trac/saga-gis/wiki/Creating%20Python%20scripts

Please don't hesitate to ask for more info or to update the wiki if you encounter bugs.


Here is a script that does the first part of the import: from las to a saga point cloud. You should proceed with converting this point cloud to a shape, perhaps after applying another module.

It uses the module io_shapes_las.dll, like explained in the wiki you can find info on the parameters by running (on windows this would be saga_cmd io_shapes_las 1).

johan@xxx$ saga_cmd libio_shapes_las 1
_____________________________________________
  #####   ##   #####    ##
 ###     ###  ##       ###
  ###   # ## ##  #### # ##
   ### ##### ##    # #####
##### #   ##  ##### #   ##
_____________________________________________

library path:   /usr/lib/saga
library name:   libio_shapes_las
module name :   Import LAS Files
author      :   O. Conrad, V. Wichmann (c) 2009
_____________________________________________
go...
Usage: 1 [-POINTS <str>] [-FILE <str>] [-T] [-i] [-a] [-r] [-c] [-u] [-n] [-R] [-G] [-B] [-e] [-d] [-p] [-C] [-VALID] [-RGB_RANGE <str>]
  -POINTS:<str>     Point Cloud
    Data Object (optional output)
  -FILE:<str>       Input File
    File path
  -T                gps-time
    Boolean
  -i                intensity
    Boolean
  -a                scan angle
    Boolean
  -r                number of the return
    Boolean
  -c                classification
    Boolean
  -u                user data
    Boolean
  -n                number of returns of given pulse
    Boolean
  -R                red channel color
    Boolean
  -G                green channel color
    Boolean
  -B                blue channel color
    Boolean
  -e                edge of flight line flag
    Boolean
  -d                direction of scan flag
    Boolean
  -p                point source ID
    Boolean
  -C                rgb color
    Boolean
  -VALID            Check Point Validity
    Boolean
  -RGB_RANGE:<str>  R,G,B value range
    Choice
    Available Choices:
    [0] 16 bit
    [1] 8 bit

I've had no means to check it because python on this computer is broken, but you get the basic concepts:

import saga_api, sys, os

##########################################
def las2spc(fLas,fSpc):
# fLas: python string containing the filename
# returns a pointer to a pointcloud
    print saga_api.SAGA_API_Get_Version()

    Library = saga_api.CSG_Module_Library()
    if not Library.Create(saga_api.CSG_String('/usr/lib/saga/libio_shapes_las.so')):
        print 'unable to load SAGA module library'
        return 0

    Module = Library.Get_Module('Import LAS Files')

    pSpc  = saga_api.SG_Create_PointCloud()
    Parms = Module.Get_Parameters() # default parameter list
    Parms('FILE').Set_Value(saga_api.CSG_String(fLas)) #filename in saga string format
    Parms('POINTS').Set_Value(pSpc) #pointcloud data object as required by the module (see above)
    Parms('T').Set_Value(0) #This are booleans, either on or off
    Parms('i').Set_Value(0) #T eg is the time, see above
    Parms('a').Set_Value(0)
    Parms('r').Set_Value(0)
    Parms('c').Set_Value(0)
    Parms('u').Set_Value(0)
    Parms('n').Set_Value(0)
    Parms('R').Set_Value(0)
    Parms('G').Set_Value(0)
    Parms('B').Set_Value(0)
    Parms('e').Set_Value(0)
    Parms('d').Set_Value(0)
    Parms('p').Set_Value(0)
    Parms('C').Set_Value(0)
    Parms('VALID').Set_Value(0)
    Parms('RGB_RANGE').Set_Value(0)

    if not Module.Execute():
        print 'module execution failed'
        return 0
    print 'module successfully executed'


    pSpc.Save(saga_api.CSG_String(fSpc))

    return pSpc 


if __name__ == '__main__':
    if len( sys.argv ) != 2:
        print 'Usage: las2shape.py <in: shape> <out: pointcloud>'
    else:
        fLas   = sys.argv[1]
        fSpc   = sys.argv[2]
        if os.path.split(fLas)[0] == '':
            fLas   = './' + fLas
        if os.path.split(fSpc)[0] == '':
            fSpc   = './' + fSpc

        pSpc = las2shape(fLas,fSpc)


#or instead of saving: use pSpc as on of the inputs to your next module; the wiki gives an example.
share|improve this answer
    
AOI = Area of Interest –  mike Sep 13 '11 at 17:59
    
That's cool ... so I don't need to call a bat file, as it seems we can run SAGA direct from python. The code isn't clear though on what tells it to use the memory rather than the HD. I have 16GB RAM. The LAS is 1km x 1km tile at 8pts/m. Looking at shp2xyz(fshp, fxyz) in the def. which line calls the SAGA command? is it if shp.Create(saga_api.CSG_String(fshp)) == 0: Are there other examples with documentation? Thanks again... –  GeorgeC Sep 14 '11 at 13:40
    
I'll try to add a more documented script later. Unless you explicitely save a file to disk it will be in memory. Eg in the first sample 01_... p('GRID') is a grid in memory. You could pass it to a different module if you want. –  johanvdw Sep 14 '11 at 14:02
    
Thanks mate. I can't download the saga_api from saga-project.org/download/python -any idea where else it is? This gives a 404 error. –  GeorgeC Sep 15 '11 at 0:46
    
That's a different project with the same name: you should look to sourceforge.net/projects/saga-gis/files/SAGA%20-%202.0/… and check for python 2.6 or 2.7 (whatever you use). I have no experience with using python on windows - better ask for help on the saga user forum –  johanvdw Sep 15 '11 at 6:58
show 3 more comments

I cannot see the code, as work stops me accessing it, but have you thought about working in the 'in_memory' workspace? It is an awful lot quicker, but do remember to clear it down.

share|improve this answer
    
Thanks that's great...i was trying to find out how to do this though and couldn't. The system also crashes on 1km x 1km LAS conversions to shp eventhough I have 1tb of space on the drives. –  GeorgeC Sep 14 '11 at 13:30
    
I can't figure out how to use the in_memory workspace to call the in memory spc file into the SAGA point cloud to shape utility. The code text is at sourceforge.net/projects/saga-gis/forums/forum/790705/topic/… and I've updated the question. –  GeorgeC Sep 15 '11 at 5:50
    
whenever you process, process to your in_memory e.g. set your out feature class to 'in_memory\\nameofoutputfclass'. remember though, that in_memory' data ceases to exist when the application stops/closes - help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… –  Hairy Sep 15 '11 at 7:53
    
-isn't this only for arcgis? I am calling a SAGA command in 'bat' file from an arcgis python script. –  GeorgeC Sep 16 '11 at 0:52
    
I tried it and it can't save it to this space from the bat. [[C:\SAGA\saga_cmd io_shapes_las "Import LAS Files" -POINTS=in_memory\\%1.spc -FILE=%3\%1_1K_class.las -T -i -a -r -c -u -n -R -G -B -e -d -p -C -VALID -RGB_RANGE=0 PAUSE c:\SAGA\saga_cmd pointcloud_tools "Point Cloud to Shapes" -POINTS=in_memory\\%1.spc -SHAPES=%2\%1_fc.shp ]] –  GeorgeC Sep 16 '11 at 1:30
show 5 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.