Latest Contributions by ChrisZemp

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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ... View more

Hi there, seem to be having trouble with the composite check(s). What check do I need to use to identify where point features are not intersecting one of two feature classes?   For example:   I have water line fittings as a point feature class.   I want to make sure all fittings are placed on (e.g. intersecting) EITHER our lateral lines feature class OR our main lines feature class.   I've tried using a composite check to create two geometry on geometry checks:    1.) Identify fittings not intersecting main lines.    2.) Identify fittings not intersecting lateral lines.   It is clear that the fitting I am testing is connected to the water lines in our geometric network, as the fitting is functioning as a node (e.g. I can click on the fitting in edit mode and click and drag it, causing the water lines it is connected with to drag with it). Yet this fitting keeps getting flagged in the reviewer table.   This isn't working. Is there something I am missing? How do I make my composite check into an OR condition? Thanks for any help or advice. ... View more