Geographic Information Systems Asked by Anna Riling on June 18, 2021
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_workspacepurple.tif")
# read in our clippers
# clipper = "D:Google DriveSchoolColor Anomaly DetectionGISBackgroundscolor_anomalycolor_anomaly.gdbortho_clip4"
bg_clipper = "D:Google DriveSchoolColor Anomaly DetectionGISBackgroundscolor_anomalycolor_anomaly.gdbbackground_clip"
targets = "D:Google DriveSchoolColor Anomaly DetectionGISBackgroundscolor_anomalycolor_anomaly.gdbtarget_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_workspaceresults_purple.csv", "w") as csv:
csv.write("index,sd,plus_minus,tp,fn,bg_anomaly,bg_not_anomalyn")
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")
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.
Answered by Sam Montoia on June 18, 2021
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.
Answered by Anna Riling on June 18, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP