I cannot get a Con expression to work. Or, if it is working, I cannot get the raster produced by the Con expression to save.
I have two rasters and I want the positive difference between the two, and if any difference is negative, for that to be set to zero.
I read on the documentation that for a spatial analysis to be performed, the rasters must be cast as raster objects, which I have done by using the arcpy.Raster(raster) tool.
import arcpy
import os
from arcpy import env
from arcpy.sa import *
Geo_output = arcpy.GetParameterAsText(0)
Geo_name = arcpy.GetParameterAsText(1)
SR = arcpy.GetParameterAsText(2)
InLAS = arcpy.GetParameterAsText(3)
pnt_spacing = arcpy.GetParameterAsText(4)
OutFC_1 = arcpy.GetParameterAsText(5)
class_1 = arcpy.GetParameterAsText(6)
return_1 = arcpy.GetParameterAsText(7)
OutFC_2 = arcpy.GetParameterAsText(8)
class_2 = arcpy.GetParameterAsText(9)
return_2 = arcpy.GetParameterAsText(10)
clip = arcpy.GetParameterAsText(11)
OutFC = [OutFC_1, OutFC_2]
classif = [class_1, class_2]
returns = [return_1, return_2]
FD_name = Geo_name + "_Feature_Dataset"
FD_addr = os.path.join(Geo_output,Geo_name + ".gdb") # Points to inside gdb
FC_addr = os.path.join(FD_addr, FD_name) # Points to inside FD
arcpy.CheckOutExtension('Spatial')
try:
arcpy.AddMessage("Creating Geodatabase: "+ Geo_name)
arcpy.CreateFileGDB_management(Geo_output, Geo_name)
arcpy.AddMessage("Creating Feature Dataset: "+ FD_name)
arcpy.CreateFeatureDataset_management(FD_addr, FD_name, SR)
for fc,cls,ret in zip(OutFC,classif,returns):
arcpy.AddMessage("Creating Multipoint Feature Class: "+ fc)
result= arcpy.LASToMultipoint_3d(InLAS,os.path.join(FC_addr,fc),pnt_spacing,
cls,ret,"#",SR,"las","1","NO_RECURSION")
name = fc + "_Terrain"
arcpy.AddMessage("Creating Terrain: "+ name)
arcpy.CreateTerrain_3d(FC_addr,name,pnt_spacing,"50000","#",
"ZTOLERANCE","ZMIN","NONE","1")
terrain = os.path.join(FC_addr,name)
arcpy.AddMessage("Adding Pyramids to: "+ name)
pyramids = ["0.5 1000"]
arcpy.AddTerrainPyramidLevel_3d(terrain,"",pyramids)
embed = fc + "_embed"
feature_class = "{} SHAPE masspoints 1 0 0.5 true false "\
"{} <None> false".format(result,embed)
arcpy.AddMessage("Adding Feature class(es) to: "+ name)
arcpy.AddFeatureClassToTerrain_3d(terrain, feature_class) # Insert Clip
arcpy.AddMessage("Building: "+name)
arcpy.BuildTerrain_3d(terrain,"NO_UPDATE_EXTENT")
raster_name = fc +"_RAS"
raster = os.path.join(FD_addr, raster_name)
bitType = "FLOAT"
method = "NATURAL_NEIGHBORS"
sampling = "CELLSIZE 10"
pyrLvl = 0
arcpy.AddMessage("Terrain to Raster: " + raster)
arcpy.TerrainToRaster_3d(terrain, raster, bitType,
method, sampling, pyrLvl)
arcpy.AddMessage("Creating Elevation Raster")
env.workspace = FD_addr
elev = os.join.path(FD_addr,"elevation")
list = arcpy.ListRasters("*", "ALL")
size = len(list) #### Make some exception here if != 2
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
**elevation = arcpy.sa.Con(raster1-raster0 >= 0, raster1-raster0,0)
elevation.save("elevation")**
except:
print arcpy.GetMessages(2)
I have triple checked to make sure that the raster list actually has files in it (it does!), and I have other rasters being saved to the FD_addr address, so I know that's not the problem either.
I'm new to python and arcgis. I'm not sure what I'm doing wrong. I've tried quite a few different ways, and none seems to work.
Help please!
Con
is not the best solution here: take the maximum of 0 and DSM-DEM instead. Using that alternative (or any alternative) will at least tell you whether the problem is withcon
or with something else.