3

I have this statement working in arcpy but I want this to implement in gdal, especially with gdal_calc, if possible.

Rastercalc = Con(((C-B) > 2) & (B > 100), A, Con(((A-D)>15) & ((A-E)>15),0,A))

I already tried writing this in gdal_calc as below but the output is not the same

gdal_calc.py --format GTiff --type Byte -A %s --A_band %s -B %s -C %s -D %s --D_band %s -E %s --E_band %s --outfile %s --calc="A*((C-B) > 2)+ (A*((A-E) < 30)* ((A-D) < 30))* ((C-B) <= 2)" --overwrite

I even tried this using numpy.where() with no luck

outcalc = numpy.where(numpy.logical_and((C - B) > 2, B > 100), A, numpy.where(numpy.logical_and((A - D) > 30, (A - E) > 30), 0, A))

Would somebody knows if this possible in gdal_calc (this is preferred way) but any other option in gdal would be great including numpy.where()?

I am pasting the part of my script and sampled results. The input data is actually a multiband (14 bands) dataset.

script:

--arcpy implementation----

            A = Raster(infile + os.sep + 'Layer_'+str(b))
            D = Raster(infile + os.sep + 'Layer_'+str(b-1))
            E = Raster(infile + os.sep + 'Layer_'+str(b+1))
            B = Raster(tempMed)
            C = Raster(tempMean)
            Rastercalc = Con(((C-B) > 2) & (B > 100), A, Con(((A-D)>15) & ((A-E)>15),0,A))
            Rastercalc.save(out_band)

            ##------ gdal_calc implementation------

            runCalc     = 'gdal_calc.py --format GTiff --type Byte -A %s --A_band %s -B %s -C %s -D %s --D_band %s \
                        -E %s --E_band %s --outfile %s --calc="(((C.astype(int16)-B > 2) * (B.astype(int16) > 100)) *\
                        ((A.astype(int16)-D <= 15) * (A.astype(int16)-E <= 15))) * A" --overwrite'\
                        % (infile, b, tempMed,tempMean,infile, b-1,infile, b+1, out_band)

spectral profiles of original layer, layer produced in arcpy, and gdal.

original

enter image description here

arcpy

enter image description here

gdal

enter image description here

Improve this question

2 Answers 2

Reset to default
5

Con(((C-B) > 2) & (B > 100), A, Con(((A-D)>15) & ((A-E)>15),0,A))

Trying to interpret the arcpy Con statement:

Where (C-B > 2) AND (B > 100) : A
Otherwise, where (A-D > 15) AND (A-E > 15) : 0
Elsewhere : A

To me this looks like it could be reduced to the following:

Where (C-B <= 2) AND (B <= 100) AND (A-D > 15) AND (A-E > 15) : 0
Elsewhere : A

Or inversely:

Where (C-B > 2) AND (B > 100) AND (A-D <= 15) AND (A-E <= 15) : A
Elsewhere : 0

This could be accomplished with the outcalc:

"((C-B > 2) * (B > 100) * (A-D <= 15) * (A-E <= 15)) * A"

The idea is each comparison yields 1 or 0 for True or False. The multiplication train acts as a progressive masking operation which only yields 1 where all conditionals are met. This value of 1 is multiplied by A to yield A and 0 times A equals 0.

Improve this answer
5
  • Thank you Logan for prompt answer with explanation but this is still not producing valid result. I even tried separating the condition in two blocks as below but still got wrong output. (((C-B > 2) * (B > 100)) * ((A-D <= 30) * (A-E <= 30))) * A
    – Devdai
    Mar 17, 2017 at 15:29
  • @Devdai What is wrong about the output? If these are Byte (unsigned 8-bit) images there is the potential for integer underflow. C-B "should" return -3 but returns 253 because negative integers can't be represented. If that seems like the reason, then you could cast the arrays in memory to a signed 16-bit for the conditional checks and recast the result to uint8. It would look something like "(((C.astype(int16)-B.astype(int16) > 2) .... *A).astype(uint8)" Watch your parentheses and cast the final result to uint8. This is within numpy namespace so don't use numpy.uint8.
    – Logan Byers
    Mar 17, 2017 at 16:06
  • Please see below for detail response. Arc produced the valid result (removing spike is the intention) and gdal produced entirely different output.
    – Devdai
    Mar 17, 2017 at 18:14
  • @Devdai I don't understand why the output image has 14 bands. It should have 1 band. It should also be a GTiff, from what your code says. You are showing a .img file. Try to delete the output file and run again. In some cases GDAL will modify existing rasters rather than write a completely new file. Also, please move your response into an edit in the original question and delete the answer you provided.
    – Logan Byers
    Mar 17, 2017 at 19:35
  • oh! my apology for the incomplete script/information. The equation is implemented to a band at a time (loop through) and later stacked using gdal_merge to match input properties. I verified the output of individual band with the final 14 band image and the result is matching.
    – Devdai
    Mar 17, 2017 at 20:20
1

I finally solved the equation. Thanks @Logan for the help. I had to break the equation to make it work.

runCalc1 =

'gdal_calc.py --format GTiff --type Byte -A %s --A_band %s -B %s -C %s -D %s --D_band %s -E %s --E_band %s --outfile %s --calc="((C-B <= 2) * (A-D >= 15) * (A-E >= 15)) * 15" --overwrite' % (infile, b, tempMed,tempMean,infile, b-1,infile, b+1, tempB1b)

runCalc2 =

'gdal_calc.py --format GTiff --type Byte -A %s --A_band %s -B %s -C %s -D %s --outfile %s --calc="((C-B > 2) * (B > 100) * (D != 15)) * A" --overwrite' % (infile, b, tempMed,tempMean, tempB1b, out_band)

Improve this answer

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.