1

So I'm coding in arcpy and trying to process all rasters in a geodatabase and output to a new geodatabase.

The steps go Minus(constant - raster), 1 divided by output, times 100. Then Reclassify results. Reclassify tool did not work. it would output a raster identical to the input.

Tried as a series of Con and Plus statements. Seemed to work. However I am running into an issue where a Con statement will ignore part of the expression for the raster and give the opposite of the desired results. This occurs when run from the native toolbox as well as from python. It has also done this whether using arcpy.gp.Con_sa or arcpy.sa.Con

ex. arcpy.gp.Con_sa(inraster,"1",savepath, 1,0,""""VALUE" > 0 AND "VALUE" <= 25"""). The input raster is a 32-bit floating point raster.

In some cases this will work correctly if the raster has values that exist between 0 and 25. However, if the raster has values that aren't between 0 and 25, it won't work correctly. The function will return a raster that has the true value for the entire raster when it should have a false value for the entire raster. It only seems to do that if the input raster values are relatively close to the expression values. I'm not really sure what to do.

Reminder input is a geodatabase containing 32-bit floating point rasters. "c" is a number input

#IMPORT
import arcpy
import os
from arcpy.sa import *
#-------------------------------------------------------------------------------------------
#ENVIRONMENT SETTINGS
#Prevents outputs from being added to map
arcpy.env.addOutputsToMap = False

#Allows intermediate layer to be overwriten
arcpy.env.overwriteOutput = True
#-------------------------------------------------------------------------------------------
#PARAMETERS
#Sets workspace as input parameter and as variable
arcpy.env.workspace = z = arcpy.GetParameterAsText(0)

#sets output workspace
destPath = arcpy.GetParameterAsText(1)

#Variables to do with naming and processing. 
a = arcpy.GetParameterAsText(2)
b = arcpy.GetParameterAsText(3)
c = arcpy.GetParameter(4)

#-------------------------------------------------------------------------------------------
#OUTPUT SAVING
#splits destination path, and gets geodatabse name and path in separate variables
x = destPath.split("\\")
for destGDB in x:
    print destGDB

destPath = destPath.strip(destGDB)

#puts the split variables back together in a new variable
finalGDB = str(destPath) + str("\\") + str(destGDB)
#-------------------------------------------------------------------------------------------
#GEODATABASE DESTINATION
#if destination geodatabase exists, do nothing, if it doesn't, make it.
if arcpy.Exists(finalGDB):
    finalGDB = finalGDB
else:
    arcpy.CreateFileGDB_management(destPath, destGDB, "CURRENT")
#-------------------------------------------------------------------------------------------
#MAIN CODE FOR LOOP

#lists rasters in input workspace
Listrasters = arcpy.ListRasters()

#Loops through rasters in input Geodatabase
for raster in Listrasters:

    #gets the layer's name in a variable
    rastname = str(raster)

    #sets final layer name
    finalname = str("Quart_") + a + str("_") + b + str("_") + rastname

    #sets the full path of the final layer
    finalpath = finalGDB + str("//") + finalname

    arcpy.AddMessage("Minusing...")
    tempras = Minus(c,raster)

    arcpy.AddMessage("Dividing...")
    tempras1 = Divide(1, tempras)

    arcpy.AddMessage("Multiplying...")
    tempras2 = Times(tempras1, 100)
    #|------------------------------------------------------------------------------------
    #Conditional statements to reclassify input

    arcpy.AddMessage("Con1...")
    workpath1 = z + str("//") + str("con1")
    con1 = arcpy.gp.Con_sa(tempras2,"1", workpath1, "0", """"VALUE" >= 0 AND "VALUE" <= 25""")

    arcpy.AddMessage("Con2...")
    workpath2 = z + str("//") + str("con2")
    con2 = arcpy.gp.Con_sa(tempras2,"2", workpath2, "0", """"VALUE" > 25 AND "VALUE" <= 75""")

    arcpy.AddMessage("Con3...")
    workpath3 = z + str("//") + str("con3")
    con3 = arcpy.gp.Con_sa(tempras2,"3", workpath3, "0", """"VALUE" > 75 AND "VALUE" < 90""")

    arcpy.AddMessage("Con4...")
    workpath4 = z + str("//") + str("con4")
    con4 = arcpy.gp.Con_sa(tempras2,"4", workpath4, "0", """"VALUE" >= 90""")

    arcpy.AddMessage("Con5...")
    workpath5 = z + str("//") + str("con5")
    con5 = arcpy.gp.Con_sa(tempras2,"5", workpath5, "0", """"VALUE" < 0""")
    #|-----------------------------------------------------------------------------------
    #Addition statments to add con rasters for full reclass

    arcpy.AddMessage("Add1...")
    add1 = arcpy.sa.Plus(con1,con2)

    arcpy.AddMessage("Add2...")
    add2 = arcpy.sa.Plus(con3,con4)

    arcpy.AddMessage("Add3...")
    add3 = arcpy.sa.Plus(con5,add1)

    arcpy.AddMessage("Add4...")
    outras = arcpy.sa.Plus(add2,add3)

    arcpy.AddMessage("Saving...")
    outras.save(finalpath)
    #|------------------------------------------------------------------------------------
    lastmsg = str("SUCCESS! You have processed ") + rastname + str("!")
    arcpy.AddMessage(lastmsg)

Can anyone advise on what's going wrong? I suspect that it may be a version issue, or an issue with the input data. But at this point I have no idea.

This issue also appears when the con tool is run on its own from its native toolbox.

Improve this question
2
  • 1
    I don't understand your methodology.. surely you could use Reclassify from arcpy.sa (resources.arcgis.com/en/help/main/10.2/index.html#//…) instead of a bunch of Con statements then adding them. Note that Raster objects implement mathematical operators so you could say tempras = (1 / ( c - raster ) ) * 100 to do all your calculations in one line. When you run each tool from the toolbox do you get the intended results? If not it's your method that's likely the problem. Note that c is likely a string and should be float(c)
    – Michael Stimson
    Feb 16, 2017 at 22:12
  • 1
    Start your script with "from arcpy.sa import *".. The rest will be easy walk. E.g. Con() because you'll be able to use temp rasters and simplified syntax. Look at the help for any of SA functions
    – FelixIP
    Feb 16, 2017 at 22:16

1 Answer 1

Reset to default
1

Your loop could be simplified thus:

import sys
import os

for raster in Listrasters:
    r         = arcpy.Raster(raster)
    rastname  = str(raster)                          #gets the layer's name in a variable
    finalname = "Quart_{}_{}_{}".format(a,b,rastname)#sets final layer name
    finalpath = os.path.join(finalGDB ,finalname)    #sets the full path of the final layer
    arcpy.AddMessage("Calculating...")
    tempras   = ( 1 / (float(c) - r) ) * 100
    arcpy.AddMessage("Properties...")
    arcpy.CalculateStatistics_management(tempras,1,1)
    MinValR = arcpy.GetRasterProperties_management(tempras,"MINIMUM")
    MaxValR = arcpy.GetRasterProperties_management(tempras,"MAXIMUM")
    MinVal = float(MinValR.getOutput(0))
    MaxVal = float(MaxValR.getOutput(0))

    arcpy.AddMessage("Reclassify...")
    reclass   = arcpy.sa.Reclassify(tempras,"value",arcpy.sa.RemapRange([
        ['%f'%MinVal,0,5],
        [0,25,1],
        [25,75,2],
        [75,90,3],
        [90,'%f'%MaxVal,4]]))
    arcpy.AddMessage("Saving...")
    reclass.save(finalpath)

As I said in my comment Raster Objects support operators, it's called Map Algebra so your whole calculation could be done in a single line. Your cascading Con and Plus statements seem to emulate a Reclassify with a RemapRange.. I have used the statistical minimum and maximum of the raster to end the ranges upper and lower as these are the limits of what values can be expected.

Give that a try and see if it helps.

Edit

Well, there's a trap for young players... it seems that it's the scientific notation in the RemapRange that it doesn't like:

>>> rng = arcpy.sa.RemapRange([[numpy.float32(MinVal),0.0,5],[0.0,25.0,1],[25.0,75.0,2],[75.0,90.0,3],[90.0,MaxVal,4]])
>>> print rng
-3.40282e+38 0.0 5; 0.0 25.0 1; 25.0 75.0 2; 75.0 90.0 3; 90.0 9.11999988556 4

but if it's specified like this:

>>> rng = arcpy.sa.RemapRange([['%f'%MinVal,0.0,5],[0.0,25.0,1],[25.0,75.0,2],[75.0,90.0,3],[90.0,'%f'%MaxVal,4]])
>>> print rng
-340282346638529010927431635313163763712.000000 0.0 5; 0.0 25.0 1; 25.0 75.0 2; 75.0 90.0 3; 90.0 9.120000 4

it works.. Thanks python for your duck types coming back to bite.

Improve this answer
10
  • When I try running your code I get TypeError: unsupported operand type(s) for -: 'float' and 'unicode' in the Map Algebra expression. Also, "c" is actually Double.
    – Lists
    Feb 16, 2017 at 23:07
  • Thank you for the quick help btw. Tried it with Minus and Divide, worked, but the reclassify threw an error. ExecuteError: ERROR 000622: Failed to execute (Reclassify). Parameters are not valid. ERROR 000628: Cannot set input into parameter remap.
    – Lists
    Feb 16, 2017 at 23:15
  • ah, I see, Unfortunately I did not test it... I'll fix it and get back to you.
    – Michael Stimson
    Feb 16, 2017 at 23:16
  • It seems that the sys.float_info isn't valid in the tool see gis.stackexchange.com/questions/153771/…. When replaced with a real value like -29999 it works just fine but that leaves scope for values outside the range. It would be nice if the values were all encompassing... I'll see what I can come up with. The first problem is the raster is actually a string and it can't find it, I have changed that to make it a Raster.
    – Michael Stimson
    Feb 16, 2017 at 23:45
  • tried it with real values. Seems to be working now. Did a stress test on a geodatabase. Works the majority of the time. but doesn't seem to process every raster. Didn't work on 6 of about 160. it seems to skip reclassifying rasters that have negative values. Will try again with your updated code. Thanks again.
    – Lists
    Feb 17, 2017 at 0:13

Your Answer

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

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