0

I'm trying to do something that is conceptually very straightforward. I'm pulling a center point from a raster (geotiff, etc) and then doing a CONTAINS operation with that point against a UTM zone grid. The intent is to obtain some field values for the polygon feature the raster center point falls within. Then, I can use those values to initiate a Project Raster method using the correct UTM zone. I have already done this with GDAL/QGIS but need to try and make it work in the Esri world.

Here is my code:

import arcpy

def getRasterCenter(raster):
    xmin = arcpy.GetRasterProperties_management(raster, 'LEFT')
    xmin = arcpy.GetRasterProperties_management(raster, 'RIGHT')
    xmin = arcpy.GetRasterProperties_management(raster, 'BOTTOM')
    xmin = arcpy.GetRasterProperties_management(raster, 'TOP')

    center_x = (float(xmax.getOutput(0)) + float(xmin.getOutput(0))) / 2
    center_y = (float(ymax.getOutput(0)) + float(ymin.getOutput(0))) / 2
    center_coords = center_x, center_y

    return center_coords

#stash intermediate data in a scratch db
tempData = arcpy.CreateScratchName(workspace=arcpy.env.scratchGDB)

#write json to file since arcpy doesn't read json string
utm_grid_geojson = '{ validated Esri json content here }'
file_json = r'C:\Temp\utm_zones.json'
with open(file_json, 'w') as outfile:
    outfile.write(utm_grid_geojson)

# set up the raster
terrain_file = r'C:\Temp\raster.tif'
terrain_file_center = getRasterCenter(terrain_file)
x = float(terrain_file_center[0])
y = float(terrain_file_center[1])
point = arcpy.Point(x, y)
terrain_center = arcpy.PointGeometry(point)

# sanity check the point
print(terrain_center.WKT)

# now get the utm grid again and do the join
utm_grid = arcpy.JSONToFeatures_conversion(file_json, tempData)

zone = arcpy.SpatialJoin_analysis(utm_grid, terrain_center, tempData, "#", "KEEP_COMMON", "#", "CONTAINS")

# check if valid result
result_desc = arcpy.Describe(zone)
print('result type: {}'.format(result_desc.featureType))

This code results in the following error:

arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid. ERROR 000670 Output Feature Class is same as input Target Features.

I am trying to avoid creating intermediate data sets if possible as this script is only one small part of a much larger model. Looking for a bit of hand-holding to understand how to do spatial joins via arcpy, I guess.

  • 1
    You are making things a bit more difficult for you. Just use the in_memory workspace for storing a proper feature class instead of working with individual geometries. Spatial Join won't be able to update geometry objects - it has been designed to work with feature layers hence the error message. – Alex Tereshenkov Oct 18 '18 at 20:35
  • Or optionally select layer by location then use a search cursor to read the values. – Michael Stimson Oct 19 '18 at 2:35
  • Thanks Alex, I went ahead and just created a feature class and it's working now, though I am astounded at how much harder it is to do programmatic reprojection with arcpy, compared to gdal.Warp() where I can just issue it the file, the EPSG code and let it go. – auslander Oct 19 '18 at 16:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.