Set spatial ref of a geoJSON polygon feature

[JamesCrandall]

Tag: Please help Jason Scheirer

After running down the path of geoJSON, I've finally encountered some difficulty in setting the spatial reference on a polygon feature.  In the code below, I thought I could get away with somehow setting the spatial ref on arcpy.AsShape() method but now I'm stuck on what to do!

The problem seems to be that my polygon.projectAs(outSr) is not working, I presume because the input polygon feature does not have the inSr set?

polygon = arcpy.AsShape(v)  # <-- does not have a spatial ref at this point
prjpolygon = polygon.projectAs(outSr)

Input JSON coordinates are web Mercator (Auxiliary Sphere).

Hopefully I'm just missing the obvious.

data = {}
feature_info = """[{"rings":[[[-9020934.359395614,3186437.9104202176],[-9020934.359395612,3186323.2548777903],[-9021029.90568097,3185711.7586515094],[-9021794.27596382,3185941.0697363648],[-9020934.359395614,3186437.9104202176]]]}]"""

def geo_convert(ring_string):  
    from json import loads, dumps  
     
    rings = loads(ring_string)  
    feat_coll = {'type': 'FeatureCollection',  
                 'features':[]}  
     
    for ring in rings:
        feat_coll['features'].append(  
            {'type': 'Feature',  
             'geometry': {  
                'type': 'Polygon',  
                'coordinates': [ring['rings'][0]]
            }})  
     
    return dumps(feat_coll)  

tmpArr = []
jsonFrmt = geo_convert(feature_info)
jsonData = json.loads(jsonFrmt)
inSr = arcpy.SpatialReference(3857)
outSr = arcpy.SpatialReference(26758)
for key,value in jsonData.iteritems():
    if value == 'FeatureCollection': 
        pass        
    else:
        for i in value:
            try:
                for k,v in i.iteritems():
                    if k == 'geometry':
                        
                        polygon = arcpy.AsShape(v) #cannot set a sr
                        
                        #project to output sr
                        prjpolygon = polygon.projectAs(outSr) #does not project

                    else:
                        pass                                        
            except:
                pass

[BruceHarold]

Hi

You need to make a geometry instance from the array you get out of the geometry feature made by AsShape().

The spatial reference is read only with any geometry object.

If your geometry collection is all polygons then AsShape(<geojsondictionary>,False) should make a multipart polygon.

If you have a big GeoJSON handling problem then consider the Data Interoperability extension, it can read GeoJSON input and write Esri geometry with less pain :-).

Regards

BTW, Jason S. no longer works here

[JamesCrandall]

You need to make a geometry instance from the array you get out of the geometry feature made by AsShape().

That's what I hoped to accomplish but I don't see how to do this with a polygon!  Do you have an example to point to?

Thanks a bunch for your input!

I can get to populate "features" array, but not sure why it won't set to a polygon

for key,value in jsonData.iteritems():
    if value == 'FeatureCollection': 
        pass        
    else:
        for i in value:
            try:
                for k,v in i.iteritems():
                    if k == 'geometry':
                        features = []
                        featureinfo = v['coordinates']
                        
                        for feature in featureinfo:
                           features.append(arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature]), inSr))

                        #project to output sr
                        prjpolygon = features.projectAs(outSr) #fails

[BruceHarold]

Hi

Here is a code snippet:

I found a multipart polygon from some polygon features I had laying around:

with arcpy.da.SearchCursor('Suburbs','shape@') as cursor:

    for row in cursor:

        gj = row[0].__geo_interface__

        if row[0].partCount > 1:

            print('Biggie')

            break

Biggie

The variable 'gj' now has GeoJSON in it.  Grab the array from geometry made from it and create a new geometry with a spatial reference

aPoly = arcpy.Polygon(arcpy.AsShape(gj,False).getPart(),arcpy.SpatialReference(2193))

Check the spatial reference:

aPoly.spatialReference.name

'NZGD_2000_New_Zealand_Transverse_Mercator'

Regards

[JamesCrandall]

Edit: Bruce -- thanks for your input and example, I'll be able to use this elsewhere I'm sure!  I should have posted the full implementation so that you didn't have to guess but I think Darren supplied the exact example I needed.

Thanks again.  I'm not understanding your example.  Where is that SearchCursor set from?  I don't have a FeatureClass, just an input JSON string that I need to take in, create the polygon feature of those coordinates, change the spatial reference and then return the coordinate pair of the centroid and the area of that feature.

[BruceHarold]

Hi, the SearchCursor was just fo rme to get hold of a feature to create GeoJSON.

I guess you're dealing with that as an input.

[JamesCrandall]

Hi Bruce,

Yes, again -- apologies for not posting my complete implementation and just offering bits and pieces.  This is all supposed to be lightweight operation and just deal with geometries minus any feature classes and cursors.

Thanks again for your input!  I learned some things!

[DarrenWiens2 · Accepted solution]

If you're going to use projectAs, it needs to be on a polygon-by-polygon basis, not the list:

>>> inSr = arcpy.SpatialReference(4326)
... outSr = arcpy.SpatialReference(3857)
... features = []
... features2 = []
... feature = [[1,2],[2,2],[2,3],[1,2]]
... polygon = arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature]), inSr)
... polygon2 = polygon.projectAs(outSr)
... features.append(polygon)
... features2.append(polygon2) 
... arcpy.CopyFeatures_management(features,r'in_memory\poly')
... arcpy.CopyFeatures_management(features2,r'in_memory\poly2')

[JamesCrandall]

Edit: the issue is an extra "[" and "]" at the start and end of the input array of coords.  Super easy to miss with all of the brackets and squiglies but your "feature" array shows the differnce!

Your example input of coords:

feature = [[1,2],[2,2],[2,3],[1,2]]  

My input looked like this:

feature = [[[-9020934.359395614, 3186437.9104202176], [-9020934.359395612, 3186323.2548777903], [-9021029.90568097, 3185711.7586515094], [-9021794.27596382, 3185941.0697363648], [-9020934.359395614, 3186437.9104202176]]]

But to remove/fix the offending brackets, I had to modify my geo_convert def that takes in the input json string and handles it.  It was an easy fix to apply to the 'coordinates' token:

def geo_convert(ring_string):  
    from json import loads, dumps  

    rings = loads(ring_string)  
    feat_coll = {'type': 'FeatureCollection',  
                 'features':[]}  

    for ring in rings:
        feat_coll['features'].append(  
            {'type': 'Feature',  
             'geometry': {  
                'type': 'Polygon',  
                'coordinates': ring['rings'][0] #<- this had an extra set of [] wrapping it
            }})  
    print feat_coll 
    return dumps(feat_coll)  
#the input string I am getting
feature_info = """[{"rings":[[[-9020934.359395614,3186437.9104202176],[-9020934.359395612,3186323.2548777903],[-9021029.90568097,3185711.7586515094],[-9021794.27596382,3185941.0697363648],[-9020934.359395614,3186437.9104202176]]]}]"""
#format the string
jsonFrmt = geo_convert(feature_info)

Darren,

I'm getting a RuntimeError "Point: Input value is not numeric" on this line:

polygon = arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature]), inSr)

Here's my input JSON:

feature_info = """[{"rings":[[[-9020934.359395614,3186437.9104202176],[-9020934.359395612,3186323.2548777903],[-9021029.90568097,3185711.7586515094],[-9021794.27596382,3185941.0697363648],[-9020934.359395614,3186437.9104202176]]]}]"""

This is what feature looks like when I print it:

[[[-9020934.359395614, 3186437.9104202176], [-9020934.359395612, 3186323.2548777903], [-9021029.90568097, 3185711.7586515094], [-9021794.27596382, 3185941.0697363648], [-9020934.359395614, 3186437.9104202176]]]

Any ideas?

features = []
features2 = []
feature = v['coordinates']
print feature

polygon = arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature]), inSr)
polygon2 = polygon.projectAs(outSr)
features.append(polygon)
features2.append(polygon2)
                       
print polygon.area
print polygon2.area

[JamesCrandall]

In case my last reply is confusing, Darren's suggested example seems to be what I needed. 

FYI for others: this is largely an issue of JSON string manipulation work and can be tough to pick out the error in the details.  Keep in mind that it's often a good idea to reduce things down to bare minimum, get it working, then apply it to the full implementation -- something easily skipped (like I tend to do!).