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?

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


# 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)
print("writing index_r")

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

# create CSV
with open("F:\\color_anomaly_workspace\\results_purple.csv", "w") as csv:
    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

  • @user2856 done! Sorry, I'm new to python. I removed some extraneous lines. – Anna Riling Aug 14 '19 at 4:16
  • @BERA It was my understanding that Extract by Mask would give you the donut "hole", not the donut "ring". Do you have a different experience? – Anna Riling Aug 14 '19 at 13:11

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.

  • Hey @Sam Montoia! My code actually has the correct indentation; I had to manually add four spaces to get this to appear as code on this site. I'll email you my latest version. – Anna Riling Aug 15 '19 at 19:25
  • Adding code doesn't seem like the most friendly of experiences here. – Sam Montoia Aug 15 '19 at 19:42
  • @AnnaRiling adding code is easy. Just paste it in, select it and click the code format button - { } – user2856 Aug 15 '19 at 21:25

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.

correctly rendered rasters without displacement

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.