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 arcpy
import pyproj as proj
import numpy as np
from geopy import distance
__author__ = 'zemp'
import arcpy, os, sys, subprocess, time, csv
from datetime import datetime
EnterpriseDB = "C:\\Users\\zemp\\AppData\\Roaming\\Esri\\Desktop10.6\\ArcCatalog\\DBO@City_Works_Test.sde"
arcpy.ImportToolbox('M:/ArcGIS/Tools/SkagitPUD.tbx', 'pud')
sqldbname = 'CityworksGIS'
sqlservername = 'arcgis'
workspace = EnterpriseDB
arcpy.env.workspace = workspace
arcpy.env.overwriteOutput = True
print "Blocking Connections..."
arcpy.AcceptConnections(EnterpriseDB, False)
print "Disconnecting Users..."
arcpy.DisconnectUser(EnterpriseDB, "ALL")
datasets = arcpy.ListDatasets("*", "Feature")
tableList = arcpy.ListTables("*", "ALL")
fcListDatasets = []
if len(datasets) > 0:
for dataset in datasets:
fcListDatasets += arcpy.ListFeatureClasses("*", "", dataset)
wIndex = [i for i, s in enumerate(datasets) if 'Water_Distribution' in s]
wIndex = int(wIndex[0])
fcListWater = arcpy.ListFeatureClasses("*", "", datasets[wIndex])
fcListRoot = arcpy.ListFeatureClasses("*", "")
if len(datasets) > 0:
fcList = fcListDatasets + fcListRoot
else:
fcList = fcListRoot
fcListWaterLength = len(fcListWater)
fcListWaterRange = range(fcListWaterLength)
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]
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]]
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")
pointList_layer = ["wCasing_layer", "wUtilityCrossing_layer", "wPipeAlignment_layer", "wRefPoint_layer", "wSystemValve_layer",
"wHydrant_layer", "wMeter_layer", "wControlValve_layer", "wFitting_layer"]
pointList_layer.sort()
nad83 = proj.Proj(init='epsg:6318')
wa_sp = proj.Proj(init='epsg:2285', preserve_units = True)
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"
memory = r"in_memory/"
tempEnd = r"_temp"
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)):
totalCount[j].append(arcpy.GetCount_management(pointList_layer[j])[0])
print("now working on",pointList_layer[j])
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
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
dist_diff = distance.geodesic(GNSS_location, GIS_location, ellipsoid='GRS-80').feet
if dist_diff >= toleranceValue:
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]) + "...
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)
edit.stopOperation()
edit.stopEditing(True)