Programmatically Clipping and Resampling an Elevation Image Service using Arcpy

[JohnGaiot]

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.

[DennisFraser]

ESRI likely does not want you to be able to make a copy.