Shift in raster in python 3.6

Question

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")

Answer 1

I wondered if your code was only applying itself to the first instance, as the indentation isn't correct for a for loop.

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

This might make the code apply to all the indexes in your index_list.

Answer 2

As usual, it was a simple thing that was missed. I had neglected to set the output coordinate system and apply a transformation. I added these lines at the beginning:

env.outputCoordinateSystem = arcpy.SpatialReference("NAD 1983 UTM Zone 13N")
env.geographicTransformations = "WGS_1984_(ITRF00)_To_NAD_1983"

The rasters then rendered properly.