1

Background: I am writing a tool to automatically project a layer (vector or raster) in ArcGIS Pro to its appropriate UTM zone.

Approach: The entire UTM zone grid is stored in a geojson string inside the script. From the string, create a temporary polygon feature class in the in_memory space, pulling only the EPSG code in as a field value. Get the centroid of the raster or vector layer. Run a spatial join between the centroid and the UTM grid. Obtain the UTM zone EPSG code. Run a Project() (if vector layer); run a ProjectRaster() (if raster layer); output the results to the project's default geodatabase.

Problem: arcpy.analysis.SpatialJoin() is not returning any rows, even though it looks like I've got correctly constructed feature classes in in_memory.

Code:

import arcpy
from arcpy.sa import *

# methods
def make_temp_utm_fc(gj_string):
    sr = arcpy.SpatialReference(4326)
    temp_fc = arcpy.management.CreateFeatureclass('in_memory', 'utmzones', 'POLYGON', spatial_reference=sr)
    fields = ['EPSG']  # list in case i want to add more attributes later
    for field in fields:
        arcpy.management.AddField(temp_fc, field, 'TEXT')

    with arcpy.da.InsertCursor(temp_fc, ['SHAPE@'] + fields) as irows:
        for feature in gj_string['features']:
            geom = arcpy.AsShape(feature['geometry'], False)
            irows.insertRow([geom] + [feature['properties']['EPSG']])

    return temp_fc

def get_extent_centroid(layer):
    desc = arcpy.Describe(layer)
    sr = desc.spatialReference
    if desc.dataType in ['FeatureLayer', 'RasterLayer']:
        extent = desc.extent
        array = arcpy.Array()
        array.add(extent.lowerLeft)
        array.add(extent.lowerRight)
        array.add(extent.upperRight)
        array.add(extent.upperLeft)
        array.add(extent.lowerLeft)
        pg = arcpy.Polygon(array, sr)
        
        return pg.centroid

# env
p = arcpy.mp.ArcGISProject('CURRENT')
arcpy.env.Workspace = p.defaultGeodatabase
arcpy.env.overwriteOutput = True

# utm zones (abbreviated for here)
utmgrid = {
    "type": "FeatureCollection",
    "name": "utmgrid",
    "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}},
    "features": [
        {"type": "Feature", "properties": {"fid": 1.0, "SWLON": "-180", "SWLAT": "-80", "HEMISPHERE": "s", "ZONE": "01", "Zone_Hemi": "01,s", "EPSG": "32501"}, "geometry": {"type": "Polygon", "coordinates": [[[#the coord ring]]]}}]}

# do the work
inlyr = arcpy.GetParameter(0)
inlyr_desc = arcpy.Describe(inlyr)

utm_fc = make_temp_utm_fc(utmgrid)

center = get_extent_centroid(inlyr)
center_row = [(center.Y, center.X)]
center_fc = arcpy.management.CreateFeatureclass("in_memory", "layercenter", "POINT", spatial_reference = inlyr_desc.spatialReference)
with arcpy.da.InsertCursor(center_fc, ['SHAPE@XY']) as cursor:
    cursor.insertRow(center_row)

# this is the part that returns nothing
sj = arcpy.analysis.SpatialJoin(center_fc, utm_fc, "in_memory/sj")

with arcpy.da.SearchCursor("in_memory/sj", "EPSG") as cursor:
    for row in cursor:
        arcpy.AddMessage(row)  # returns None

Once I have the EPSG code from the containing UTM zone polygon, it's easy to do the reprojection so I haven't included that code. Where am I going wrong doing the spatial join?

Currently I'm working exclusively with layers that carry a 4326 spatial reference, so both the centroid and UTM grid are in GCS_WGS_1984. I've added AddMessage() calls throughout for debugging and it seems like it's successfully created the UTM grid and centroid as Polygon and Point feature classes respectively.

3
  • I'm personally unfamiliar with using geojson in arcpy. I would start there and confirm if your code is actually creating a valid geometry and EPSG field, try writing it to a shapefile rather than in_memory and then check it in ArcPro? In fact write everything to a file geodatabase, get your code working then update to write to in_memory, this is often what I do as in_memory can sometimes be a black box. – Hornbydd Apr 26 at 18:56
  • So, simply by moving outputs to the default geodatabase, everything worked without issue. Which is kind of just as odd, because it means that the in_memory space is failing to hold/overwrite the feature classes for some unknown, inscrutable reason. (it also takes about triple the time for it to complete with the added I/O operations). – auslander Apr 26 at 19:27
  • I have never been able to confirm if data in in_memory have spatial indices. So something like a spatial join with probably run slower despite the data being in a faster data store... – Hornbydd Apr 26 at 21:25

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.