I have a geoJSON service I'm querying that's being supplied by GeoServer. I have a python script that ingests this feed routinely and uses arcpy.AsShape to render the geoJSON as various feature classes (I am extracting points, lines and polygons).
The problem is that the X/Ys in the service itself are backwards, even according to its own documentation. So my features are falling in all kinds of impossible locations.
Is there a fast way in python to flip the X/Y geometry in feature classes, including the vertices of line and polygon types?
I'm not aware of something like a "flip" geoprocessing tool, but you can do this using arcpy geometry objects and an Update Cursor.
Is the service public? If so, it would be helpful to share a link.
Creating and manipulating geometries objects is more expensive than creating and manipulating native Python data types. Instead of using arcpy.AsShape() and then trying to flip X/Y, why not process the incoming GeoJSON service and do the flipping before passing to arcpy.AsShape() to create the geometry?
Joshua,
Unfortunately the service is not public. However I can give you an idea of the payload. It looks like this:
{
"type":"FeatureCollection",
"totalFeatures":4000000,
"features":
[
{
"type":"Feature",
"id":(guid here)
"geometry":
{
"type":"Point",
"coordinates":[10.0000000000,11.0000000000]
}
...
and for "LineString" type, it reads
"coordinates":[[10.0000000000,11.0000000000],[12.0000000000,13.0000000000],etc....]
So in every. single. list item, I need to swap the lat/lons.
I am barely getting out of novice python and am not very savvy when it comes to this. I am of course interested in doing this in the least costly way, so doing this prior to saving it as a series of feature classes is fine; I just figured there may be something very brief built into arcpy that might do this without much glue on an extant feature class.
If this is not the case I'm looking for help either with a snippet of code or even pseudocode to help me think through the mechanics.
It should be mentioned that currently I am ingesting this content as a .json file and then writing feature classes out of it. This content represents millions of rows, so I am not sure what would be quicker - opening an existing feature class, or re-writing the text document. Of course I am totally open to being shown better ways to do this!
Before I dive into it any deeper, I want to understand how you are processing this GeoJSON information. Are you taking the JSON file and iterating over the feature collection and calling arcpy.AsShape() with each feature/geometry?
Joshua, yes, that's correct. It looks like this:
with arcpy.da.InsertCursor(pointFC, ['SHAPE@XY'] + fields) as prows:
for feature in geojson['features']:
geom = arcpy.AsShape(feature['geometry'], False)
if geom.type == 'point':
prows.insertRow([geom] + feature['properties']
using the SHAPE@ token for line/poly types after that, etc..
Not incorporating the cursor part, the following will go through the features and swap all the coordinate pairs for points, lines, etc.... and create geometries:
for feature in geojson['features']:
geom_geojson = feature['geometry']
coords = geom_geojson['coordinates']
if geom_geojson['type'] == "Point":
geom_geojson['coordinates'] = coords[::-1]
else:
geom_geojson['coordinates'] = [coord[::-1] for coord in coords]
geom = arcpy.AsShape(geom_geojson, False)
>>> a = np.arange(2*5).reshape(2,5).T
>>> a
array([[0, 5],
[1, 6],
[2, 7],
[3, 8],
[4, 9]])
>>> np.flipud(a)
array([[4, 9],
[3, 8],
[2, 7],
[1, 6],
[0, 5]])
>>> np.fliplr(a)
array([[5, 0],
[6, 1],
[7, 2],
[8, 3],
[9, 4]])
easier as well
Somehow, I find this is easier to get my head around...
>>>
>>> List = [[1,2],[3,4], [5,6], [6,7]]
>>> for l in List:
... l.reverse()
...
>>> List
[[2, 1], [4, 3], [6, 5], [7, 6]]
>>>