1

I'm very new to python and arcpy.

The task: I have a large series of rasters (r0, r1, r2, r3...) (tif; n=240) and corresponding quality rasters (qr0, qr1, qr2, qr3...) (tif) with values 0,1,2,3 (whereas 0 is good quality and all other values are bad quality). I would like to check the cell quality of a quality raster (e.g. qr1) in a if-statement. If condition is true (quality=good=0), I would like to replace that quality value by the corresponding raster (r1), if not, I would like to check if the quality is true in the previous quality raster (qr0) at the same pixel location and if true, I would like to replace the quality value with that raster value (r0). In the case that quality value is also bad, I would like to check a third case: if the following quality raster value (qr2) is true, than replace the value with the corresponding raster (r2). If the quality is bad in all three quality rasters (qr0, qr1, qr2), than assign a no data to the cell.

My initial idea: As far as I'm informed I should use a nested for loop-statements with if or elif-statements but I'm a bit confused in how assigning the speficic rasters and looping through them.

  • 1
    Did you check this link gis.stackexchange.com/questions/163368/… ? You can use Raster Calculator for this purpose. – Shiko Nov 21 '16 at 10:04
  • I already saw this link before. But my difficulty is in looping through the workspace and assigning the previous and following quality raster and raster with replace values. I could solve it exemplary for one sequence (see following expression) but I would like to loop through all files. Con("qr1.tif" == 0,"r1.tif", Con("qr0" == 0,"r0",Con("qr2.tif" == 0,"r2.tif", -999))) . Furthermore, I'm willing to learn and understand python. – xmisx Nov 21 '16 at 10:40
2

SCRIPT:

import arcpy, time, os
from arcpy import env
from arcpy.sa import *
env.overwriteoutput=True
## PARAMETERS
qualityGDB=r'...\Scratch.gdb'
dataGDB=r'D:\SCRATCH'
# GET RASTERS. MAKE EFFORTS TO ADJUST THEM IF THEIR ORDER MISMATCH !
arcpy.env.workspace = qualityGDB
qRasters = arcpy.ListRasters("*", "GRID")
arcpy.env.workspace = dataGDB
dRasters = arcpy.ListRasters("*", "GRID")

for i in range(2,len(dRasters)):
    dList,qList=[],[]
    for j in range(i,i-3,-1):
        dList.append(arcpy.Raster(dataGDB+os.sep+dRasters[j]))
        qList.append(arcpy.Raster(qualityGDB+os.sep+qRasters[j]))
        arcpy.AddMessage(j)
    pick=Con(CellStatistics(qList,"MINIMUM") == 0,LowestPosition(qList))
# BUG FOUND Pick doesn't work as expected
##    answer=Pick(pick, dList)
    arcpy.gp.Pick_sa(pick, dList, r"D:/aerials/images%s%s.tif" %(os.sep,str(i).zfill(3)))

SAMPLE:

enter image description here

| improve this answer | |
  • I had and still have some problems with my Python-Version so I couldn't test it. I followed the suggestion of Shiko and used the Raster Calculator for the time being. But thanks for your proposal. I will test it very soon! – xmisx Nov 29 '16 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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