UpdateCursor trouble... modified geometry must be a different instance

Question

Hi there. I have created a script that works as a GPS integrity script for all of our point features. If a user has inadvertently moved a point feature from it's GPS'ed location, than the script will move it back, keeping our accuracy high. In other words, the script calculates the distance between the 'SHAPE@X' and 'SHAPE@Y' fields and our high accuracy GNSS enabled fields "ESRIGNSS_LONGITUDE" and "ESRIGNSS_LATTITUDE". If the calculated distance (dist_diff) fall outside the bounds of our set tolerance level (toleranceValue) for that particular feature class, then the points are moved back to our GNSS locations via the UpdateCursor. The GNSS lat and lon points are re-projected to our local state plane coordinate system using PyProj. Then it should be as simple as replacing the SHAPE@X and SHAPE@Y fields that fall outside the toleranceValue with the re-projected lat and lon values, and then writing those features to a new feature class. However, I keep getting the following error: Traceback (most recent call last): File "", line 184, in RuntimeError: The modified geometry must be a different geometry instance from the feature's original geometry (e.g., a copy or new instance). [class name = CityworksGIS_Test.DBO.wControlValve, object id = 104] This is puzzling to me, as the first feature class in the list passes through the UpdateCursor with no problems. But then I get this error on the second iteration of the loop, with the wControlValve feature class. I feel like I could just be overlooking something... If someone could please put their eyes on my code, I'd truly appreciate it, and would gladly return the favor in the future. Thank you. Code below: # Import These import arcpy import pyproj as proj import numpy as np from geopy import distance ### ENIVRONMENT SETUP ### __author__ = 'zemp' import arcpy, os, sys, subprocess, time, csv from datetime import datetime # Variables # Enterprise Database Connection EnterpriseDB = "C:\Users\zemp\AppData\Roaming\Esri\Desktop10.6\ArcCatalog\DBO@City_Works_Test.sde" # Import Toolbox arcpy.ImportToolbox('M:/ArcGIS/Tools/SkagitPUD.tbx', 'pud') # The name of the SQL database for script logging with pyodbc sqldbname = 'CityworksGIS' # The name of the SQL Server sqlservername = 'arcgis' # Set workspace workspace = EnterpriseDB # Set the workspace environment arcpy.env.workspace = workspace arcpy.env.overwriteOutput = True # block new connections to the working and prod database. print "Blocking Connections..." arcpy.AcceptConnections(EnterpriseDB, False) # disconnect all users from CityworksGIS database. print "Disconnecting Users..." arcpy.DisconnectUser(EnterpriseDB, "ALL") # -------------------------------------DATA IMPORT------------------------------------------------------------------ # Gather all datasets datasets = arcpy.ListDatasets("*", "Feature") # Gather all tables tableList = arcpy.ListTables("*", "ALL") # Gather all feature classes fcListDatasets = [] if len(datasets) > 0: for dataset in datasets: fcListDatasets += arcpy.ListFeatureClasses("*", "", dataset) # Get wFCs from water distribution dataset wIndex = [i for i, s in enumerate(datasets) if 'Water_Distribution' in s] wIndex = int(wIndex[0]) fcListWater = arcpy.ListFeatureClasses("*", "", datasets[wIndex]) # Gather FCs from ROOT: fcListRoot = arcpy.ListFeatureClasses("*", "") # Combine Lists if len(datasets) > 0: fcList = fcListDatasets + fcListRoot else: fcList = fcListRoot # -------------------------------------FIND INDICES AND FILEPATHS OF DATA------------------------------------------ fcListWaterLength = len(fcListWater) fcListWaterRange = range(fcListWaterLength) # Find index of feature classes within feature class list: hydrantIndex = [i for i, s in enumerate(fcListWater) if 'wHydrant' in s] meterIndex = [i for i, s in enumerate(fcListWater) if 'wMeter' in s] mainLineIndex = [i for i, s in enumerate(fcListWater) if 'wMain' in s] controlValveIndex = [i for i, s in enumerate(fcListWater) if 'wControlValve' in s] fittingIndex = [i for i, s in enumerate(fcListWater) if 'wFitting' in s] systemValveIndex = [i for i in fcListWater if 'wSystemValve' in i and '_D' not in i] casingIndex = [i for i in fcList if 'CasingEndpoint_SurveyGrade' in i] utilityCrossingIndex = [i for i in fcList if 'UtilityCrossing_SurveyGrade' in i] pipeAlignmentIndex = [i for i in fcList if 'PipeAlignment_SurveyGrade' in i] refPointIndex = [i for i in fcList if 'ReferencePoint_SurveyGrade' in i] # Pull file paths wCasing = casingIndex[0] wUtilityCrossing = utilityCrossingIndex[0] wPipeAlignment = pipeAlignmentIndex[0] wRefPoint = refPointIndex[0] wSystemValve = workspace + "\" + systemValveIndex[0] wHydrant = workspace + "\" + fcListWater[hydrantIndex[0]] wMeter = workspace + "\" + fcListWater[meterIndex[0]] wControlValve = workspace + "\" + fcListWater[controlValveIndex[0]] wFitting = workspace + "\" + fcListWater[fittingIndex[0]] # Combine into list pointList = [wCasing, wUtilityCrossing, wPipeAlignment, wRefPoint, wSystemValve, wHydrant, wMeter, wControlValve, wFitting] arcpy.MakeFeatureLayer_management(wCasing, "wCasing_layer") arcpy.MakeFeatureLayer_management(wUtilityCrossing, "wUtilityCrossing_layer") arcpy.MakeFeatureLayer_management(wPipeAlignment, "wPipeAlignment_layer") arcpy.MakeFeatureLayer_management(wRefPoint, "wRefPoint_layer") arcpy.MakeFeatureLayer_management(wSystemValve, "wSystemValve_layer") arcpy.MakeFeatureLayer_management(wHydrant, "wHydrant_layer") arcpy.MakeFeatureLayer_management(wMeter, "wMeter_layer") arcpy.MakeFeatureLayer_management(wControlValve, "wControlValve_layer") arcpy.MakeFeatureLayer_management(wFitting, "wFitting_layer") # Combine into list pointList_layer = ["wCasing_layer", "wUtilityCrossing_layer", "wPipeAlignment_layer", "wRefPoint_layer", "wSystemValve_layer", "wHydrant_layer", "wMeter_layer", "wControlValve_layer", "wFitting_layer"] pointList_layer.sort() ### METHOD 1: ### nad83 = proj.Proj(init='epsg:6318') # assuming you're using NAD83 wa_sp = proj.Proj(init='epsg:2285', preserve_units = True) # use a locally appropriate projected CRS gnssX = [[] for x in range(len(pointList))] gnssY = [[] for x in range(len(pointList))] shapeX = [[] for x in range(len(pointList))] shapeY = [[] for x in range(len(pointList))] objectID = [[] for x in range(len(pointList))] distDiff = [[] for x in range(len(pointList))] totalCount = [[] for x in range(len(pointList))] outPath = [[] for x in range(len(pointList))] toleranceValue = 0 flag = "_Flagged" # Path for writing FCs to memory: memory = r"in_memory/" tempEnd = r"_temp" # Path for new Feature Dataset fd = "C:\Users\zemp\AppData\Roaming\Esri\Desktop10.6\ArcCatalog\DBO@City_Works_Test.sde\GPSIntegrityCheck" edit = arcpy.da.Editor(workspace) edit.startEditing(False, True) edit.startOperation() for j in range(len(pointList_layer)): # Get count for each point feature totalCount[j].append(arcpy.GetCount_management(pointList_layer[j])[0]) print("now working on",pointList_layer[j]) # Set tolerance Values if 'wHydrant' in str(pointList_layer[j]): toleranceValue = 2.0 featureName = "wHydrant" + flag elif 'wControlValve' in str(pointList_layer[j]): toleranceValue = 1.0 featureName = "wControlValve" + flag elif 'wSystemValve' in str(pointList_layer[j]): toleranceValue = 0.5 featureName = "wSystemValve" + flag elif 'wFitting' in str(pointList_layer[j]): toleranceValue = 0.5 featureName = 'wFitting' + flag elif 'wMeter' in str(pointList_layer[j]): toleranceValue = 1.0 featureName = "wMeter" + flag elif 'CasingEndpoint' in str(pointList_layer[j]): toleranceValue = 0 featureName = 'wCasingEndPoint' + flag elif 'UtilityCrossing' in str(pointList_layer[j]): toleranceValue = 0 featureName = 'wUtilityCrossing' + flag elif 'PipeAlignment' in str(pointList_layer[j]): toleranceValue = 0 featureName = 'wPipeAlignment' + flag elif 'ReferencePoint' in str(pointList_layer[j]): toleranceValue = 0 featureName = 'wReferencePoint' + flag else: pass # Start search cursor for the following rows print("The tolerance value is", toleranceValue, "for", pointList_layer[j]) with arcpy.da.UpdateCursor(pointList_layer[j],['SHAPE@X', 'SHAPE@Y', 'ESRIGNSS_LONGITUDE', 'ESRIGNSS_LATITUDE', 'OBJECTID']) as cursor: for row in cursor: if row[2] is not None and row[3] is not None: [lonX, latY] = proj.transform(wa_sp, nad83, row[0], row[1]) GNSS_location = row[3],row[2] GIS_location = latY, lonX # Calculate difference between GIS and GNSS location using NAD83 and GRS Elipsoid dist_diff = distance.geodesic(GNSS_location, GIS_location, ellipsoid='GRS-80').feet if dist_diff >= toleranceValue: # If points is outside tolerance zone, do this: objectID[j].append(row[4]) gnssX[j].append(lonX) gnssY[j].append(latY) shapeX[j].append(lonX) shapeY[j].append(latY) distDiff[j].append(dist_diff) row[0] = lonX row[1] = latY cursor.updateRow(row) else: pass else: pass print(pointList_layer[j], "The total count is: " + str(totalCount[j]) + "...# of records outside tolerance: " + str(len(distDiff))) for id in objectID[j]: stringID = str(id) where_clause = "OBJECTID = '" + stringID + "'" arcpy.SelectLayerByAttribute_management(pointList_layer[j], "ADD_TO_SELECTION", where_clause=where_clause) if j != 0: del totalSelected totalSelected = int(arcpy.GetCount_management(pointList_layer[j])[0]) else: totalSelected = int(arcpy.GetCount_management(pointList_layer[j])[0]) print ("New Selection Total: " + str(totalSelected) + " Selection from tolerance table: " + str(len(objectID[j])) + " Records that were not matched: " + str(len(objectID[j])-totalSelected) + " Total Records: " + str(totalCount[j][0])) outPath = os.path.join(fd, str(pointList_layer[j])+"_GPS_Flagged") arcpy.CopyFeatures_management(pointList_layer[j], outPath) # Stop Editing edit.stopOperation() edit.stopEditing(True)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

ChrisZemp · Answer

Hi Josh,

Thanks for your reply. I tried all variations of the SHAPE@ token and talked with American ESRI support. They told me it's been a known issue since 2012, and that the behavior of the SHA{E@ token is inconsistent, primarily dependent on your working GDB environment (e.g. SDE vs localized GDB, etc). They basically told me it wasn't possible.

I did some more google-fu and some more brain storming. I changed my code and figured out how to pythoncially drop and re-create our geometric network without causing any hiccups or errors with our SDE and Cityworks integration. I'm currently working on a blog post about my solution and will post it here when done, in case anyone else runs into this error in the future. Cheers

JoshuaBixby · Answer

I typically work with the whole shape token (SHAPE@), or sometimes XY (SHAPE@XY).  I suspect this may be related to tolerances and using SHAPE@X and SHAPE@Y.  I suggest changing your code to use SHAPE@ to see if replacing the point geometry with a completely new point geometry makes the error go away.