Community

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

722
2
08-05-2020 10:52 AM
ChrisZemp
by
New Contributor II
‎08-05-2020 10:52 AM

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 "<input>", line 184, in <module>
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) + "\nSelection from tolerance table: " + str(len(objectID[j]))
      + "\n Records that were not matched: " + str(len(objectID[j])-totalSelected) + "\n 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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
2 Replies
JoshuaBixby
by MVP Esteemed Contributor
MVP Esteemed Contributor
‎08-05-2020 11:57 AM

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.

0 Kudos
ChrisZemp
by
New Contributor II
‎10-13-2020 09:24 AM

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