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. I have tried using a Con(IsNull()) statement as well as ExtractByMask to create the "donut hole" where the target is in the background raster. These methods work, except that the background 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 this step 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.
Why am I seeing a shift in this one and only raster, and how do I prevent it?
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 clippers
# clipper = "D:\\Google Drive\\School\\Color Anomaly Detection\\GIS\\Backgrounds\\color_anomaly\\color_anomaly.gdb\\ortho_clip4"
bg_clipper = "D:\\Google Drive\\School\\Color Anomaly Detection\\GIS\\Backgrounds\\color_anomaly\\color_anomaly.gdb\\background_clip"
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")
# create CSV
with open("F:\\color_anomaly_workspace\\results_purple.csv", "w") as csv:
csv.write("index,sd,plus_minus,tp,fn,bg_anomaly,bg_not_anomaly\n")
print("wrote csv")
# add indices to a list and then begin loop
index_list = [r, g]
print("begin loop")
for index in index_list:
# 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.50, 1.1, 0.5)
for multiplier in multipliers:
high_threshold = mean + multiplier * sd
low_threshold = mean - multiplier * sd
# set extent
arcpy.env.extent = "MAXOF"
arcpy.env.snapRaster = index
arcpy.env.cellSize = index
# 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_target = ExtractByMask(diff_rc, targets_rast)
diff_rc_target.save(basename + multiplier_name + "_target.tif")
# convert background clip polygon to raster
arcpy.AddField_management(bg_clipper, "raster_id", "SHORT")
arcpy.CalculateField_management(bg_clipper, "raster_id", 1, "PYTHON")
bg_clipper_rast = arcpy.PolygonToRaster_conversion(bg_clipper, "raster_id", "background_clip_rast.tif")
# create background raster
diff_rc_background = ExtractByMask(diff_rc, bg_clipper_rast)
diff_rc_background.save(basename + multiplier_name + "_background.tif")
# build attribute table to get cell counts
arcpy.BuildRasterAttributeTable_management(diff_rc_target)
arcpy.BuildRasterAttributeTable_management(diff_rc_background)
arcpy.CheckInExtension("spatial")