Latest Contributions by JohnGaiot

‎07-21-2025

That does not make sense to me. Then why make an entire raster function module/library available in Python when you are not allowed to use it? This is our own data which is open without restrictions, so we should be able to process it to meet our own internal requirements.

‎07-10-2025

I am a GIS Data Analyst working for the Ontario provincial government and we have a series of public lidar elevation image services available here. I am working with an unrestricted version of one of our lidar image services in our Dev environment at a native resolution of approximately 0.5m resolution. Ultimately, I would like to programmatically resample a portion to a coarser resolution for further processing. Our organization is currently running ArcGIS Image Server 10.9.1. Specifically, I would like to clip/extract a portion of the image service data programmatically in Arcpy by watershed - a relatively small area for now at 700 sq.km. but eventually would like to go bigger. Then resample the selection to 2m resolution and work with that derivative product in further processing. I've managed to connect to the service and tried using Spatial Analyst tools like ExtractByMask but they take too long to process as expected (beyond an overnight process which we had to bail out of). Ideally, I would like to use raster functions on the service through Arcpy. I could find very little if any documentation on this, specifically relating to examples of using the "arcgis.raster.functions" module or libaries. In ArcGIS Pro, they are recognized and import fine as shown below. I create an ImageServiceLayer connection to our Dev image service, and then try to apply the clip raster function on it as a first step. import arcpy, os, sys import arcgis.raster.functions as rf from arcgis.raster import Raster from arcpy.sa import * # Check out the Spatial Analyst extension license arcpy.CheckOutExtension("Spatial") # Functions Definition Section def is_running_in_python_window(): return arcpy.GetArgumentCount() == 0 def AddMsg(InText): if is_running_in_python_window(): print(InText) else: AddMsg(InText) # Input Parameters OutputWS = arcpy.GetParameterAsText(0) if OutputWS == '#' or not OutputWS: OutputWS =r"C:\_Workspace\_Data\OHN\PredictedStreams\FlintRiver_TestRun2\NewTest\ElevationSourceTest" # provide a default value if unspecified HardCodedService = arcpy.GetParameterAsText(5) # Hardcoded service option for elevation input if HardCodedService == '#' or not HardCodedService: HardCodedService = "true" # provide a default value if unspecified Input_DTM = arcpy.GetParameterAsText(4) # Source elevation at native lidar or imagery resolution if Input_DTM == '#' or not Input_DTM: ProcessDTM = False if HardCodedService == "true": AddMsg("Note: Using hardcoded DTM service for elevation input..") Input_DTM = r"https://lrc-ndmnrf-webr.test.gov.on.ca/arcgis/rest/services/Elevation/FRI_DTM_SPL_Trial/ImageServer" # Internal provincial service for 0.5m resolution single photon lidar ProcessDTM = True else: ProcessDTM = True # Input_DTM = r"\FRI_DTM_SPL_Trial" # provide a default value if unspecified Output_Buffered_SA = arcpy.GetParameterAsText(7) if Output_Buffered_SA == '#' or not Output_Buffered_SA: Output_Buffered_SA = OutputWS + "\\VectorProcessing.gdb\\Buffered_SA" # provide a default value if unspecified Output_ResampledDTM = arcpy.GetParameterAsText(10) if Output_ResampledDTM == '#' or not Output_ResampledDTM: Output_ResampledDTM = OutputWS + "\\Layers\\Resampled_DTM.tif" # provide a default value if unspecified # Create a SpatialReference object for Ontario Lambert CSRS ontario_lambert_csrs = arcpy.SpatialReference(3162) # "NAD_1983_CSRS_Ontario_Lambert" # Set the environment's output coordinate system arcpy.env.outputCoordinateSystem = ontario_lambert_csrs # Print the coordinate system to verify AddMsg(f"Output Coordinate System: {arcpy.env.outputCoordinateSystem.name}") AddMsg("Setting study area extents to buffered study area..") # Get the current extents desc = arcpy.Describe(Output_Buffered_SA) polyextent = desc.extent arcpy.env.extent = polyextent if ProcessDTM: AddMsg("Resetting extents to nearest 2m..") # Round the extents to the nearest 2 meters xmin = round(polyextent.XMin / 2) * 2 xmax = round(polyextent.XMax / 2) * 2 ymin = round(polyextent.YMin / 2) * 2 ymax = round(polyextent.YMax / 2) * 2 # Create a new extent object with the rounded values AddMsg(f"Processing extents: xmin:{xmin} xmax:{xmax} ymin:{ymin} ymax:{ymax}") new_extent = arcpy.Extent(xmin, ymin, xmax, ymax) # Apply the new extent (this example assumes you are using it for a spatial reference or analysis) arcpy.env.extent = new_extent arcpy.env.cellAlignment = "ALIGN_WITH_PROCESSING_EXTENT" arcpy.env.resamplingMethod = "NEAREST" AddMsg("Processing DTM...") ProcessService = True # Set to True if using an image service, False if using a local raster if ProcessService: arcpy.env.workspace = OutputWS image_service_layer = arcpy.MakeImageServerLayer_management(Input_DTM, "ImageServiceLayer") in_raster = image_service_layer descR = arcpy.Describe(in_raster) inputCellSize = descR.meanCellHeight AddMsg(" > Clipping raster image service...") clip_raster = rf.clip(raster=in_raster, geometry=Output_Buffered_SA) # Clip the raster to the buffered study area AddMsg(" > Reprojecting...") rpj_raster = rf.reproject(clip_raster, ontario_lambert_csrs, 2, 2) # Reproject to Ontario Lambert CSRS AddMsg(" > Resampling to 2m...") resample_raster = rf.resample(raster=clip_raster, resampling_type="Minimum", input_cellsize=inputCellSize, output_cellsize=2) # Resample to 2m resolution AddMsg(" > Saving output...") resample_raster.save(Output_ResampledDTM) When I try using the 'rf.clip' raster function, it produces the following error. Note: Using hardcoded DTM service for elevation input.. Output Coordinate System: NAD_1983_CSRS_Ontario_MNR_Lambert Setting study area extents to buffered study area.. Resetting extents to nearest 2m.. Processing extents: xmin:856280 xmax:912580 ymin:12573830 ymax:12627000 Processing DTM... > Clipping raster image service... Traceback (most recent call last): File "<string>", line 96, in <module> File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\site-packages\arcgis\raster\functions\__init__.py", line 2027, in clip function_chain_ra = copy.deepcopy(template_dict) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 146, in deepcopy y = copier(x, memo) ^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 231, in _deepcopy_dict y[deepcopy(key, memo)] = deepcopy(value, memo) ^^^^^^^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 146, in deepcopy y = copier(x, memo) ^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 231, in _deepcopy_dict y[deepcopy(key, memo)] = deepcopy(value, memo) ^^^^^^^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 172, in deepcopy y = _reconstruct(x, memo, *rv) ^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 271, in _reconstruct state = deepcopy(state, memo) ^^^^^^^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 146, in deepcopy y = copier(x, memo) ^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 231, in _deepcopy_dict y[deepcopy(key, memo)] = deepcopy(value, memo) ^^^^^^^^^^^^^^^^^^^^^ File "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\copy.py", line 167, in deepcopy raise Error( copy.Error: un(deep)copyable object of type <class 'geoprocessing server result object'> Am I accessing the image service layer incorrectly? Is there an alternative or better way to utilize raster functions via Arcpy? Many thanks in advance.

‎03-22-2024

HOLY COW TI ! Your solution, although a roundabout one, worked pretty slick! I can't believe it's been so hard to find a SIMPLE solution.. Cheers.

‎03-22-2024

Why do we have to perform all this coding when simply importing a saved layer file (lyr or lyrx) should work programmatically? In other words, what's the point of having the "Apply Symbology From Layer" function when it's not doing what it's supposed to be doing?

‎02-16-2024

It seems like a basic function but in my experience, it's not. There is one function that comes very close and it's a quick one-step solution. It's called "Raster Domain".. However, the Downsides: it requires 3D Analyst, and more importantly, the final polygon goes through the centroid of each cell at the data extent of the raster, so it's not an exact 'footprint' as it falls a half cell short on all sides of features represented. But if you want something quick and approximate, this would be the ideal solution. For something more exact, the simplest method would be to use Raster Calculator (which requires a Spatial or Image Analyst license). Set the "NoData" Environment variable to something other than Zero. If you have a floating point or non-discrete raster like a DEM, you will first need to convert your raster to integer. Second, multiply the raster by Zero, so you get a single value for the entire raster where there were valid values. Third, convert the Zero raster to a polygon. Now you have your extents based on features in a raster.

‎12-17-2023

I was wondering if this issue was ever resolved. I am getting the same problem. The original raster tiles can be displayed as unique integer values but when combined into a mosaic, the unique value option is no longer available even though the bit depth is the same as before (unsigned 8 bit). Is there a work around?

‎01-20-2021

Thanks Jon for the update. Looking forward to the fix. As a workaround, I keep adding the disable/enable combination for each network created, which seems to be doing the trick. I found a bigger issue regarding a performance hit in Pro of the DA Update Cursor when processing data in loops compared to ArcMap. I will have to report this in the appropriate forum when I have a chance as it was enough to halt my development and turn to other pressing matters in the interim.

‎11-18-2020

UPDATE #3: Even stranger finding... I tried to remove the enable topology functions and it logically came up with an error that topology was not enabled when trying to set flow direction. So basically, I found after creating each network, I had to first disable and then enable topology before setting flow direction and performing a trace successfully. Very weird behaviour, but I seem to be getting results now, which are not exactly the same as I was getting in ArcMap but at least I can troubleshoot my trace results now. Below attached is a code snippet showing 1 of the 3 network/traces being performed and the weird logic that made it run successfully.

‎11-16-2020

Thanks for fixing my problem. Unfortunately, with it resolved, I can no longer access this other 'new' account so I can't flag the issue as resolved... ;(

‎11-09-2020

UPDATE: I tried saving the results of the first trace to a dataset in a new file geodatabase, created a new network (wc_net) and still got the same 'table already exists' error when I try to enable the network. Nothing else exists in the geodatabase because I just created it. So I have no clue why it's doing this...

‎11-03-2020

Thanks Jonathan, this sounds very promising and actually very similar to the way the older geometric network tracing functioned. For now, I am trying to utilize the AGGREGATED_GEOMETRY option by exploding the results of my first trace and using that in a second network. But when I try to enable it, I get a 'table already exists' error, which I suspect is related to the first network I created/enabled. Details are in a more recent post I made in this forum a few days ago. When I get a chance, I may try a test of creating the second network in a brand new geodatabase to see if that works any better?

‎10-30-2020

Hello, I am trying to port over an 10.3.1 script which used the old geometric network functionality in Python scripting. I am running ArcGIS Pro 2.6.2. In my script, I am creating multiple trace networks to set up a cleaned watercourse network for further hydrology and DEM processing (burning and flow directions etc). It creates the first network (wc_clean_net) to find connected segments (isolated segments are dropped) and enables the network no problem. It also successfully creates a second network (wc_net) in the same dataset (called "Layers") using an upstream trace this time from a particular pour point at the edge of the study area. But when I try to enable the network, I get the following error: 160326: The table already exists. Traceback (most recent call last): File "C:\_Workspace\Toolboxes\MNRF_TerrainToolsLocal\Pro\TerrainScripts\Step A.2. Prepare Water_new.py", line 756, in <module> arcpy.EnableNetworkTopology_tn(wc_net, "", "ENABLE_TOPO") File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\tn.py", line 174, in EnableNetworkTopology raise e File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\tn.py", line 171, in EnableNetworkTopology retval = convertArcObjectToPythonObject(gp.EnableNetworkTopology_tn(*gp_fixargs((in_trace_network, max_number_of_errors, only_generate_errors), True))) File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 511, in <lambda> return lambda *args: val(*gp_fixargs(args, True)) arcgisscripting.ExecuteError: ERROR 160326: The table already exists. Failed to execute (EnableNetworkTopology). Failed to execute (CreateGeometricNetwork).‍‍‍‍‍‍‍‍‍‍‍‍ I can enable it manually by running the Toolbox command when I load the network into a map but I can't do it through python. Below is a sample of the script. First I try to clean up previous runs of the script deleting the networks before creating them again to help try isolating the problem. I am also checking for "TN_" tables (related to the trace network) that might not be cleaned up properly from previous bailed runs. As a last ditch effort, I tried copying all input feature classes and datasets to a new geodatabase to try running on a 'clean' geodatabase. No luck. Line 58 is where it bails. # Network clean-up from previous runs arcpy.AddMessage("Checking for Networks from previous run...") if arcpy.Exists(wc_clean_net): arcpy.AddMessage(" >> Deleting 'wc_clean_net'...") arcpy.Delete_management(wc_clean_net) if arcpy.Exists(wc_net): arcpy.AddMessage(" >> Deleting 'wc_net'...") arcpy.Delete_management(wc_net) arcpy.Delete_management(WC_Clean) if arcpy.Exists(wc_select_net): arcpy.AddMessage(" >> Deleting 'wc_select_net'...") arcpy.Delete_management(wc_select_net) arcpy.Delete_management(WC_Select) #arcpy.env.workspace = Output_Vector_GDB_Workspace arcpy.AddMessage("Workspace = " + arcpy.env.workspace) CurrentLayer = "TN_*" arcpy.AddMessage("Checking for '" + CurrentLayer + "' orphan Network Tables...") GDB_Tables = arcpy.ListTables(CurrentLayer, "All") tbl_ctr = 0 for tbl in GDB_Tables: tbl_ctr += 1 arcpy.Delete_management(tbl) arcpy.AddMessage(" >> " + str(tbl_ctr) + " '" + CurrentLayer + "' tables deleted.") arcpy.AddMessage("Performing Connected Trace to create WC Clean layer...") WC_Name = os.path.basename(WC_Dissolve2) arcpy.AddMessage(" >> Dissolved WC Name = '" + WC_Name + "'") arcpy.AddMessage(" >> Creating Trace Network...") arcpy.CreateTraceNetwork_tn(Layers, "wc_clean_net", "", [[WC_Name, "SIMPLE_EDGE"]]) arcpy.AddMessage(" >> Enabling Network...") arcpy.EnableNetworkTopology_tn(wc_clean_net, "", "ENABLE_TOPO") arcpy.AddMessage(" >> Snapping buffer and pour points to network...") arcpy.Snap_edit(wc_bf_and_pour_pts, [[WC_Dissolve2, "EDGE", "0.1 Meters"]]) arcpy.AddMessage(" >> Tracing network...") # Process: Trace Geometric Network ##arcpy.TraceGeometricNetwork_management(wc_clean_net, WC_Selection3, wc_bf_and_pour_pts, "FIND_CONNECTED", "", "", "", "", "", "NO_TRACE_ENDS", "", "", "", "AS_IS", "", "", "", "AS_IS") arcpy.Trace_tn(wc_clean_net, "CONNECTED", wc_bf_and_pour_pts, "", "", "", "", "", "", "", "", "", "EDGES_ONLY", "", "", "AGGREGATED_GEOMETRY", "NEW_SELECTION", "", "Connected_Network_1", "", wc_traceout) arcpy.MultipartToSinglepart_management(wc_traceout, WC_Clean) wc_clean_lyr = "wc_clean_lyr" arcpy.MakeFeatureLayer_management(WC_Clean, wc_clean_lyr) arcpy.AddMessage(" >> Retrieving connected network results...") arcpy.AddMessage(" >> Layer name is: {0}. Count is: {1}.".format(os.path.basename(WC_Clean), str(arcpy.GetCount_management(wc_clean_lyr)))) arcpy.AddMessage("Performing Upstream Trace to create initial WC Select layer...") try: WC_Name = os.path.basename(WC_Clean) arcpy.AddMessage(" >> Cleaned WC Name = '" + WC_Name + "'") arcpy.AddMessage(" >> Creating Trace Network...") arcpy.CreateTraceNetwork_tn(Layers, "wc_net", "", [[WC_Name, "SIMPLE_EDGE"]]) except: arcpy.AddMessage(" >> Network exists... bypassing creation") arcpy.AddMessage(" >> Enabling Network...") arcpy.EnableNetworkTopology_tn(wc_net, "", "ENABLE_TOPO") arcpy.AddMessage(" >> Snapping pour points to network...") arcpy.Snap_edit(Flags_or_Sources_for_Tracing, [[WC_Clean, "EDGE", "0.1 Meters"]]) arcpy.AddMessage(" >> Tracing network...") arcpy.Trace_tn(wc_net, "UPSTREAM", Flags_or_Sources_for_Tracing, Barriers_for_Tracing, "", "", "", "", "", "", "", "", "EDGES_ONLY", "", "", "AGGREGATED_GEOMETRY", "NEW_SELECTION", "", "Upstream_Network_2", "", wc_traceout2) if Barriers_for_Tracing != "": arcpy.MultipartToSinglepart_management(wc_traceout2, Output_Vector_GDB_Workspace + "\\" + wc_inside) wc_inside_lyr = "wc_inside_lyr" arcpy.MakeFeatureLayer_management(Output_Vector_GDB_Workspace + "\\" + wc_inside, wc_inside_lyr) arcpy.AddMessage(" >> Retrieving upstream 'inside' network results...") arcpy.AddMessage(" >> Layer name is: {0}. Count is: {1}.".format(wc_inside, str(arcpy.GetCount_management(wc_inside_lyr)))) else: arcpy.MultipartToSinglepart_management(wc_traceout2, WC_Select) wc_select_lyr = "wc_select_lyr" arcpy.MakeFeatureLayer_management(WC_Select, wc_select_lyr) arcpy.AddMessage(" >> Retrieving upstream 'select' network results...") arcpy.AddMessage(" >> Layer name is: {0}. Count is: {1}.".format(os.path.basename(WC_Select), str(arcpy.GetCount_management(wc_select_lyr)))) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ All I can think of is that the Enable function is trying to use the same system table name for the second network that was created for the first network. Am I missing something or is this a bug? Thanks, John

‎10-27-2020

Similar to Ben Coakley's issue about "Function Outputs", there is practically no documentation that I can find about how to utilize the results of a trace (SELECTION option) using arcpy scripting. Basically, I want to save the selection to a layer so I can do further tracing or analysis on the derived selection. There should be an option similar to AGGREGATED_GEOMETRY which saves NON-AGGREGATED results. I don't see a rationale for aggregating geometry anyway as it produces just a single mult-part feature which seems useless to work with. Will exploding the aggregated geometry always maintain flow direction of the individual features? Even so, if there were important attributes on the original network, they would now be lost in the aggregated dissolve, which means you can't do any attribute queries or further selections on the aggregated output. By comparison, the old TraceGeometricNetwork function allowed you to work with the Group Layer output and save the results to a separate feature class with the attributes intact. I would like to see this functionality in the new Trace function (or at least the ability to save the selection with original geometry and attributes as optional output).

‎06-25-2020

I would argue the printed image is not correct. I would say it's more related to the rendering capability of Pro/ArcMap when trying to handle projections and maybe scale. Logic might dictate that images would appear at their best in their native projection, but oddly that's not the case. Using the ESRI World Imagery in its native WGS 84 Web Mercator projection, the image appears blurry exactly like Theodore's example above. The blurry phenomenon happens at various scales. For example, it's blurry at 1:8000 (and coarser scales) and seems best when viewed at 1:6000. However, when I change the map projection to NAD 83 CSRS UTM Zone 17N, it appears clear within a larger range of scales (between 1:3000 and 1:10000). Note I was able to replicate this in traditional ArcMap Desktop 10.3.1 (I know.. ancient!) as well. So this rendering issue does not seem to be platform specific. Perhaps it could be the graphics card?? My specs are the following although I personally think it's still a pretty good card: https://www.nvidia.com/content/dam/en-zz/Solutions/design-visualization/documents/75509_DS_NV_Quadro_K620_US_NV_HR.pdf

‎07-18-2019

Sorry Jacob.. I may have not responded directly to you.. see below or my comment earlier today on this thread.

Online Status	Offline
Date Last Visited	a week ago

My Ideas

Latest Contributions by JohnGaiot

Re: Programmatically Clipping and Resampling an Elevation Image Service using Arcpy

Programmatically Clipping and Resampling an Elevation Image Service using Arcpy

Re: Updating Data Source in ArcGIS PRO project with arcpy

Re: Change symbology of an existing layer through arcpy

Re: How to Get Polygon of Raster Extent in ArcGIS Pro

Re: How to symbolize unique values for single-band mosaiked rasters?

Re: Trace Network Enabling Problem

Re: Trace Network Enabling Problem

Re: This New Platform Update Messed Up My Account

Re: Trace Network Enabling Problem

Re: How do you work with Trace Results?

Trace Network Enabling Problem

How do you work with Trace Results?

Re: Blurry basemap in Pro layout view

Re: Calling Arc Hydro Tools in Python in ArcGIS Pro?

Re: Change symbology of an existing layer through ...

Programmatically Clipping and Resampling an Elevat...

Re: Blurry basemap in Pro layout view

Re: This New Platform Update Messed Up My Account

Re: Why can you no longer display disabled Python ...