-1

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!

10
  • You need to define a workspace for that method of .save to work.
    – Aaron
    Jul 1, 2013 at 4:08
  • @Aaron I defined it as env.workspace = FD_addr towards the end of the code, but before the .save. Is this not the correct way to do it? Is there a different/better way to save?
    – SAC_Phy
    Jul 1, 2013 at 4:17
  • Sorry, I missed that. Try saving by a direct pathname (e.g. elevation.save("C:/temp/elevation.img") ). Also try simplifying your Con function to make sure there is nothing wrong with your parameters (e.g. elevation = arcpy.sa.Con(raster1, 10, 1, "VALUE > 1"). If the results save, there is a problem with your Con function
    – Aaron
    Jul 1, 2013 at 4:33
  • @Aaron I tried both of those. Neither saving to my desktop nor simplifying my Con expression worked. I'm starting to think this is less of a "not saving" problem and more of a "not being created" problem.
    – SAC_Phy
    Jul 1, 2013 at 5:49
  • 2
    (1) My point is that by isolating the problem in your code you will make it much easier for people to read it and help you. (2) 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 with con or with something else.
    – whuber
    Jul 1, 2013 at 21:51

1 Answer 1

1

I often have trouble saving raster objects when using python scripts, which I do quite a bit. I have developed a habit of defining the workspace immediately before the raster.save, as I have found this helps debug. Have a look at the code below and see if it helps.

arcpy.AddMessage("Creating Elevation Raster")
arcpy.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 = Raster(list[0])
raster1 = Raster(list[1])
difference = raster1 - raster0
elevation = Con(difference >= 0, difference, 0)
arcpy.env.workspace = FD_addr
elevation.save("elevation")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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