Select to view content in your preferred language

# calculating coordinates from a table, non spatial process.

967
8
07-18-2023 01:02 PM
by
Occasional Contributor III

I have a table to process that with coordinates for records expressed in two different ways - one set are GCS, the other are State Plane coordinates. I'd like to simply calculate new coordinates into a new field so that there is one common set of coordinates for all records. I can do this through a multistep process involving converting recordsets of each coordinate type to spatial layers, projecting one of them, merging them back together and calculating the coordinate values. These data ultimately will become spatial, but in an early stage of processing, it would be nice to have a simpler (aka fewer steps) process. I suspect i could go find a formula for this somewhere along the way and implement that in Python, but this seems like it could be more of a pain that the multistep process i outlined above.

Tags (2)
1 Solution

Accepted Solutions
by Anonymous User
Not applicable

You could step through the features and use pyproj to transform. Not sure how your data is set up or anything really, but add two fields and use an update cursor.

import arcpy
import pyproj

fc = r''

# Convert coordinate system and add XY Coordinates in WGS 1984
with arcpy.da.UpdateCursor(fc, ['GCS_x', 'GCS_y', 'sp_x', 'sp_y', 'converted_lon', 'converted_lat']) as uCur:
for row in uCur:
# Already in GCS?, but if needed use the pyproj step below
if row[0]:
row[4] = row[0]
row[5] = row[1]
elif row[2]:
# Stateplane to WGS84
stplne = pyproj.Proj(init='epsg:32133')
wgs84 = pyproj.Proj(init='epsg:4326')
row[4], row[5] = pyproj.transform(stplne, wgs84, row[2], row[3])

uCur.updateRow(row)

8 Replies
MVP Regular Contributor

The only issue is that it doesn't do it in place.

by
Occasional Contributor III

Thanks, I fussed with this approach a bit and  still puts me in the position of having to create intermediate tables and layers. Appreciate the suggestion! I found another post that more or less outlines how to go about doing these via formula but it looked pretty complex and lots of room for error.

by
Occasional Contributor III

Do you want to perform the projection in arcpy (gdb's) or arcgis (SEDF)?  If arcpy, you could write to a memory gdb instead of gdb directly.  This might save some of your concern about extra tables.

Alternatively, you could use the arcgis API for python and a spatially enabled dataframe (SEDF).   Two possible suggestions:
1) One SEDF.  Use pandas apply to pass the geometry project method on each geometry object.  This should work even if you have mixed starting spatial references.  You can read the in_sr with the SpatialReference method into a side column.  You might have to create geometry objects first with similar pandas apply.  Docs at https://developers.arcgis.com/python/api-reference/arcgis.geometry.html#project.

2) Multiple SEDFs, one SEDF per spatial reference group of objects.  Use the GeoSeriesAccessor to set the sr of the entire SEDF shape column then use the project_as method to your desired sr.  Finally, concatenate the SEDFs.  Docs at https://developers.arcgis.com/python/api-reference/arcgis.features.toc.html#geoseriesaccessor.

HTH
Tyler

by
Occasional Contributor III

Thanks Tyler! This is work on an enterprise GDB before the layer gets published so the SEDF probably isn't my solution here, but that's an interesting approach. I need to learn more about the API for Python so maybe this will give me an excuse to poke around at it.

ct

by Anonymous User
Not applicable

You could step through the features and use pyproj to transform. Not sure how your data is set up or anything really, but add two fields and use an update cursor.

import arcpy
import pyproj

fc = r''

# Convert coordinate system and add XY Coordinates in WGS 1984
with arcpy.da.UpdateCursor(fc, ['GCS_x', 'GCS_y', 'sp_x', 'sp_y', 'converted_lon', 'converted_lat']) as uCur:
for row in uCur:
# Already in GCS?, but if needed use the pyproj step below
if row[0]:
row[4] = row[0]
row[5] = row[1]
elif row[2]:
# Stateplane to WGS84
stplne = pyproj.Proj(init='epsg:32133')
wgs84 = pyproj.Proj(init='epsg:4326')
row[4], row[5] = pyproj.transform(stplne, wgs84, row[2], row[3])

uCur.updateRow(row)

by
Occasional Contributor III

ah, I should have known there is a library for this, I'll check that out. I've already written the code to create two spatial layers, perform a projection on one then bring them back together, but using this library may be useful in simplifying the process.

New Contributor III

Hi,

I recommend strongly pyproj.

by
Occasional Contributor III

Thanks to both Tomasz and JeffK - I've implemented the pyproj module to handle this in my table and it works perfectly well. Ultimately, I think this will help streamline some downstream processes and just makes everything more straightforward. One thing I am noticing is when i run pyproj I get some warnings that I think are saying the transform method is deprecated, so the module available to import from the ArcGIs Pro package manager may be out of date. I have to read those warnings a bit closer and chase that down.

Edit: something to note here if you use this approach - the deprecated method illustrated by JeffK works, but the method shown by Tomasz is significantly faster. Working with 10K records using an insert cursor the first method took 5-6 minutes to run, the second method took under a minute, and if you read the 'gotchas' section in the pyproj docs, that holds up given some timing reported there.