print "Importing modules..."
import os
import arcpy, sys, traceback, exceptions
try:
#STEP 1: SET OVERWRITE OUTPUT TO TRUE
arcpy.env.overwriteOutput = True
print "Overwrite Output = True"
#STEP 2: CHECK OUT EXTENSION
#STEP 3: ACTIVATE THE SPATIAL ANALYST EXTENSION
print "The code is importing modules"
if arcpy.CheckExtension("3D"):
if arcpy.CheckOutExtension("3D"):
#STEP 4: SET THE WORKSPACE
arcpy.env.workspace = r"C:\Python_ArcGIS\LAS"
thePath = r"C:\Python_ArcGIS\LAS"
LAS = r"C:\Python_ArcGIS\LAS\LAS.gdb"
# Create the file gdb that will store the feature dataset
if arcpy.Exists("LAS.gdb"):
arcpy.Delete_management("LAS.gdb")
arcpy.management.CreateFileGDB(thePath, "LAS.gdb")
gdb = '{0}/{1}'.format(thePath, LAS)
#Create Feature Dataset
# CreateFeatureDataset_management (out_dataset_path, out_name, {spatial_reference})
if arcpy.Exists("LAS_FD"):
arcpy.Delete_management("LAS_FD")
arcpy.management.CreateFeatureDataset(LAS, "LAS_FD", "C:\Python_ArcGIS\SpatialReference.gdb\NAD_1983_UTM_Zone_15N")
fd = '{0}/{1}'.format(LAS, "LAS_FD")
LAS_FD = r"C:\Python_ArcGIS\LAS\LAS.gdb\LAS_FD"
# Create a LOOP for to reach all the files in the folder
for root, dirs, files in os.walk(thePath):
print "**********************"
print "ROOT: " + root
for aFile in files:
print root + "\\" + aFile
#Creating multipoint features from LAS
# LASToMultipoint_3d (input, out_feature_class, average_point_spacing, {class_code}, {return}, {attribute}, {input_coordinate_system}, {file_suffix}, {z_factor}, {folder_recursion})
extension1 = os.path.splitext(aFile)[0]
print extension1
LAS_File = extension1 + ".las"
if aFile == LAS_File:
MultiName = "DTM_MTP_" + aFile
arcpy.LASToMultipoint_3d (aFile, MultiName, 0.7, 2)
print MultiName
extension = os.path.splitext(MultiName)[0]
print extension
Multi_FC = extension + ".shp"
#Move Shapefiles to Geodatabase
#FeatureClassToGeodatabase_conversion (Input_Features, Output_Geodatabase)
arcpy.FeatureClassToGeodatabase_conversion (Multi_FC, LAS_FD)
#Create Terrain Dataset for each file
else:
print "No More Files to Process: Move to Terrain"
print "Creating Terrain"
#Create Terrain
# CreateTerrain_3d (in_feature_dataset, out_terrain_name, average_point_spacing, {max_overview_size}, {config_keyword}, {pyramid_type}, {windowsize_method}, {secondary_thinning_method}, {secondary_thinning_threshold})
arcpy.CreateTerrain_3d (LAS_FD, "DTM_Terrain", 0.7)
print "Adding terrain pyramid levels"
# Add Terrain Pyramid Level
# AddTerrainPyramidLevel_3d (in_terrain, pyramid_level_definition, {pyramid_type})
pyramid_level = ["1.5 1000", "3 2500", "6 5000", "12 10000", "24 20000", "48 50000"]
in_terrain = LAS_FD + "\\" + "DTM_Terrain"
print in_terrain
arcpy.AddTerrainPyramidLevel_3d (in_terrain, 'WINDOWSIZE', pyramid_level)
print "Creating a list of feature classes"
arcpy.env.workspace = LAS_FD
fcList = arcpy.ListFeatureClasses()
params = "fcList shape mass points <None> 0 48 <none> false <none> <none>"
print "Adding feature classes to the terrain"
#Add feature class to terrain
# AddFeatureClassToTerrain_3d (in_terrain, in_features)
# [[in_features, height_field, SF_type, group, min_resolution, max_resolution, overview, embed, embed_name, embed_fields, anchor],...]
arcpy.AddFeatureClassToTerrain_3d (in_terrain, params)
## print "Building Terrain"
## #Build Terrain
## # BuildTerrain_3d (in_terrain, {update_extent})
## arcpy.BuildTerrain_3d (in_terrain, "NO_UPDATE_EXTENT")
##
## #terrain to raster
## # TerrainToRaster_3d (in_terrain, out_raster, {data_type}, {method}, {sample_distance}, {pyramid_level_resolution})
## arcpy.TerrainToRaster_3d ("Terrain", in_terrain, "INT", "NATURAL_NEIGHBORS")
## print "All files are done"
else:
print "The 3D is not check out."
arcpy.CheckInExtension("3D")
else:
print "The 3D extension is not present."
except:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n " + str(sys.exc_info()[1])
msgs = "ARCPY ERRORS:\n" + arcpy.GetMessages(2) + "\n"
arcpy.AddError(msgs)
arcpy.AddError(pymsg)
print msgs
print pymsg
arcpy.AddMessage(arcpy.GetMessages(1))
print arcpy.GetMessages(1)
