Shift in raster using Con(IsNull()) in python 3.6

by Anna Riling   Last Updated August 14, 2019 03:22 AM - source

I have a python script for a color anomaly detection study. I am attempting to create two rasters each with three classes: one raster representing a "target" and the other representing "background", or all pixels except those within the target extent. The "target" raster was created by converting a polygon to raster. I am using a Con(IsNull()) statement, which works great, creating the "donut hole" where the target is in the background raster--except that the raster is shifted a few feet to the northwest. The target raster is in the right spot, however. I have created other rasters in the script prior to these two and all are in the correct location geographically. I set the extent, snapRaster, and cellSize to an original "spectral index" raster. Entire code is below. My questions are why am I seeing a shift in this one and only raster and how do I prevent it?

Resulting background raster is created correctly, but shifted to the northwest a few feet

import arcpy
from arcpy import env
from arcpy.sa import *
import numpy as np
from arcpy import os

# To allow overwriting the outputs change the overwrite option to true.
env.overwriteOutput = True

arcpy.CheckOutExtension("spatial")

# define main directory
main_dir = "F:\\color_anomaly_workspace\\"
env.workspace = main_dir

# read in main ortho dataset
ortho = Raster("F:\\color_anomaly_workspace\\purple.tif")

# read in our clipper
clipper = "D:\\Google Drive\\School\\Color Anomaly Detection\\GIS\\Backgrounds\\color_anomaly\\color_anomaly.gdb\\ortho_clip4"
targets = "D:\\Google Drive\\School\\Color Anomaly Detection\\GIS\\Backgrounds\\color_anomaly\\color_anomaly.gdb\\target_footprint"

# clip the mosaic to the clipper extent
ortho_clip = arcpy.Clip_management(ortho, "241896.71980000008 4137828.0887 241956.74800000008 4137967.415", "purple_clip.tif")

# create indices
basename = arcpy.Describe(ortho_clip).catalogPath
b1 = Raster(basename + "\\Band_1")
b2 = Raster(basename + "\\Band_2")
b3 = Raster(basename + "\\Band_3")

r = Float(b1)
r.save("index_r")
print("writing index_r")

g = Float(b2)
g.save("index_g")
print("writing index_g")

b = Float(b3)
b.save("index_b")
print("writing index_b")

rn = r / (r + g + b)
rn.save("index_rn")
print("writing index_rn")

gn = g / (r + g + b)
gn.save("index_gn")
print("writing index_gn")

bn = b / (r + g + b)
bn.save("index_bn")
print("writing index_bn")

m = (r + b) / 2
m.save("index_m")
print("writing index_m")

c = (b + g) / 2
c.save("index_c")
print("writing index_c")

y = (g + r) / 2
y.save("index_y")
print("writing index_y")

mn = m / (m + c + y)
mn.save("index_mn")
print("writing index_mn")

cn = c / (m + c + y)
cn.save("index_cn")
print("writing index_cn")

yn = y / (m + c + y)
yn.save("index_yn")
print("writing index_yn")

ndrg = (r - g) / (r + g)
ndrg.save("index_ndrg")
print("writing index_ndrg")

ndgb = (g - b) / (g + b)
ndgb.save("index_ndgb")
print("writing index_ndgb")

ndbr = (b - r) / (b + r)
ndbr.save("index_ndbr")
print("writing index_ndbr")

ndmc = (m - c) / (m + c)
ndmc.save("index_ndmc")
print("writing index_ndmc")

ndcy = (c - y) / (c + y)
ndcy.save("index_ndcy")
print("writing index_ndcy")

ndym = (y - m) / (y + m)
ndym.save("index_ndym")
print("writing index_ndym")

# create CSV
with open("F:\\color_anomaly_workspace\\results_purple.csv", "w") as csv:
csv.write("index,sd,plus_minus,tp,fp,tn,fn\n")

# add indices to a list and then begin loop
index_list = [r, g, b, rn, gn, bn, m, c, y, mn, cn, yn, ndrg, ndgb, ndbr, ndmc, ndcy, ndym]

for index in index_list:
# set extent
env.extent = index
env.snapRaster = index
env.cellSize = index

# calculate focal mean of spectral index
index_focal = FocalStatistics(index, NbrRectangle(600, 600, "CELL"), "MEAN")
basename = arcpy.Describe(index).basename
index_focal.save(basename + "_focal")

# perform image difference
diff = index_focal - index
diff.save(basename + "_diff")

# get mean and SD
mean = float(arcpy.GetRasterProperties_management(diff, "MEAN").getOutput(0))
sd = float(arcpy.GetRasterProperties_management(diff, "STD").getOutput(0))
min1 = float(arcpy.GetRasterProperties_management(diff, "MINIMUM").getOutput(0))
max1 = float(arcpy.GetRasterProperties_management(diff, "MAXIMUM").getOutput(0))

# begin reclassification loop
multipliers = np.arange(0.5, 5.1, 0.5)
for multiplier in multipliers:
    high_threshold = mean + multiplier * sd
    low_threshold = mean - multiplier * sd

    # reclassify raster
    rc = [[min1, low_threshold, 1],
          [low_threshold, high_threshold, 2],
          [high_threshold, max1, 3]]

    diff_rc = Reclassify(diff, "Value", RemapRange(rc))
    diff_rc.save(basename+ "_diff_rc.tif")

    # select out just purple target polygons
    targets_purple = arcpy.Select_analysis(targets, "target_purple", "color = 'purple'")
    arcpy.AddField_management(targets_purple, "raster_id", "SHORT")
    arcpy.CalculateField_management(targets_purple, "raster_id", 1, "PYTHON")
    targets_rast = arcpy.PolygonToRaster_conversion(targets_purple, "raster_id", "targets_purple_rast.tif")

    # clip raster to target extent
    multiplier_name = str(round(int(multiplier * 100),0)).rjust(3,"0")
    diff_rc_background = Con(IsNull(targets_rast), diff_rc)
    diff_rc_background.save(basename + multiplier_name + "_background.tif")
    diff_rc_target = Con(targets_rast, diff_rc)
    diff_rc_target.save(basename + multiplier_name + "_target.tif")

    # build attribute table to get cell counts
    arcpy.BuildRasterAttributeTable_management(diff_rc_target)
    arcpy.BuildRasterAttributeTable_management(diff_rc_background)

    # write results to CSV
    with open("F:\\color_anomaly_workspace\\results_purple.csv", "a") as csv:
        csv.write(basename + "," + str(multiplier) + "," + "plus" + ",")

        # process target for positive SD
        diff_rc_target_path = os.path.join(env.workspace, (basename + multiplier_name + "_target.tif"))
        sc = arcpy.da.SearchCursor(diff_rc_target_path, ["VALUE", "COUNT"])
        for row in sc:
            if row[0] == 3:
                tp = row[1]
            elif row[0] == 2:
                fn = row[1]
        sc = None

        # process background for positive SD
        diff_rc_background_path = os.path.join(env.workspace, (basename + multiplier_name + "_background.tif"))
        sc = arcpy.da.SearchCursor(diff_rc_background_path, ["VALUE", "COUNT"])
        for row in sc:
            if row[0] == 3:
                fp = row[1]
            elif row[0] == 2:
                tn = row[1]

        sc = None

        csv.write(str(tp) + "," + str(fp) + "," + str(tn) + "," + str(fn) + "\n")

        csv.write(basename + "," + str(multiplier) + "," + "minus" + ",")

        # process target for negative SD
        diff_rc_target_path = os.path.join(env.workspace, (basename + multiplier_name + "_target.tif"))
        sc = arcpy.da.SearchCursor(diff_rc_target_path, ["VALUE", "COUNT"])
        for row in sc:
            if row[0] == 1:
                tp = row[1]
            elif row[0] == 2:
                fn = row[1]
        sc = None

        # process background for negative SD
        diff_rc_background_path = os.path.join(env.workspace, (basename + multiplier_name + "_background.tif"))
        sc = arcpy.da.SearchCursor(diff_rc_background_path, ["VALUE", "COUNT"])
        for row in sc:
            if row[0] == 1:
                fp = row[1]
            elif row[0] == 2:
                tn = row[1]

        sc = None

        csv.write(str(tp) + "," + str(fp) + "," + str(tn) + "," + str(fn) + "\n")

arcpy.CheckInExtension("spatial")


Related Questions



Calculate percentile value of raster

Updated February 25, 2019 15:22 PM


set extent for jpeg programmatically

Updated September 11, 2016 08:09 AM