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 (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.
- What do all the sys.arg[1-4] values mean? Is this designed to run as an arcgis tool?
- "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.
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.
