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.

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!

share|improve this question
    
You need to define a workspace for that method of .save to work. –  Aaron Jul 1 '13 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 '13 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 '13 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 '13 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 '13 at 21:51

1 Answer 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")
share|improve this answer

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.