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
arcpy
gdal