Reprojecting tabular data without creating feature class?

533
10
01-29-2019 03:44 PM
Arne_Gelfert
Regular Contributor

Okay, I realize I may draw the scorn from the community for asking the way I'm asking but ...Is there an easy way (easier than what I'm doing) to read Lat/Lon values stored in a non-spatial SQL table, reproject them and write the new values back to same or other table? Since this is for a non-GIS audience, I don't have any use for creating an intermediate feature class. I found arcpy.AddXY_management. But I believe, as most other scenarios I could come up with, requires a feature class. Basically, I'm curious if there is another way to do this without creating a feature class.

This is what I tried:

import pyodbc

sql_conn = pyodbc.connect('''bla bla bla''') 
query = '''SELECT [LATITUDE],[LONGITUDE]
           FROM MyTable'''

cursor = sql_conn.cursor()
cursor.execute(query)

ptDataRows = []

for row in cursor.fetchall():
    #Turn cursor row into list
    ptDataRows.append([x for x in row])

I was bringing in additional attributes from my source table which I'm leaving out in this example. But a longer list would be the only difference. Next I'm working with the Lat/Lon columns to create point geometries and then 

import arcpy

src_sr = arcpy.SpatialReference(4267) 
dest_sr = arcpy.SpatialReference(4326)

for i,row in enumerate(ptDataRows):

    pt = arcpy.Point(row[0],row[1])
    ptGeom = arcpy.PointGeometry(pt, src_sr)
 
    newpoint = ptGeom.projectAs(dest_sr,'NAD_1927_To_WGS_1984_4')
    newcoord = [round(newpoint.centroid.X,7),round(newpoint.centroid.Y,7)]
 
    pt_attrs = row.extend(newcoord)

Now, I have a new list of lists, each containing the data for a new table row that I can write back to a target table. -- I also looked at arcpy.project_management. But besides the hilarious name that will make any PMP laugh, I couldn't get it to work. I think it requires a feature class again.

Since this works as it, I have to mention that I'm needing to append additional spatial attributes such as UTM zones. That's why I'm wondering if this is the "best" approach.

0 Kudos
10 Replies
DanPatterson_Retired
MVP Esteemed Contributor

If you are working with point geometries, your approach is about as easy as it gets, unless you want to start talking about adding event layers to maps and the whole Define Projection thing.  It may be easier on all involved if you roll out a toolbox with a tool that does the process behind the scenes

JoeBorgione
MVP Esteemed Contributor

Waaay back in the ArcInfo days, there was a tool that would project a file as long as you could specify what projection the x y fields were in.  I've never understood why that was never ported to ArcGIS and/or Pro...

That should just about do it....
0 Kudos
Arne_Gelfert
Regular Contributor

that's interesting. Probably way before my time. I just think it would be handy tool to have - one to just turn XY (old) into XY (new) without creating a feature class. I liked Dan's suggestions of creating a tool. Maybe a geoprocessing web tool? POST XY(old) and get XY(new) returned.

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

I agree with others, there isn't any other approach that is more straightforward or simple.  You could go code golfing and condense what you have down more, but I don't see that changing performance or results.  Since you are already using ArcPy, I suggest using ArcPy cursors to interact with the data instead of pyodbc.

Arne_Gelfert
Regular Contributor

hadn't even occurred to me to try and use arcpy coursors with anything other than feature classes. Spent some time away from GIS and was using pyodbc successfully doing other things. So the familiar module came in handy. Will have to try cursors next time.

0 Kudos
RandyBurton
MVP Regular Contributor

Cursors or just printing (don't get x, y, lon, lat mixed up):

import arcpy

latlon = arcpy.SpatialReference(4326) # wgs_1984
webmer = arcpy.SpatialReference(3857) # web mercator
nad  = arcpy.SpatialReference(4267) # nad_1927

def convertXY(x,y,inSR,outSR):
    try:
        ptGeo = arcpy.PointGeometry(arcpy.Point(x,y),inSR).projectAs(outSR)
        return (ptGeo.centroid.X, ptGeo.centroid.Y)
    except:
        print "\nERROR - cannot convert point\n  x: {}  y: {}  inSR: {}  outSR: {}".format(x,y,inSR.name,outSR.name)
        return (0.0,0.0)

tbl = r"C:\Path\to\file.gdb\XYtable"
flds = [ 'x_4326','y_4326','x_3857','y_3857','x_4267','y_4267' ]

# using an update cursor, lon and lat fields are populated, calculate others
with arcpy.da.UpdateCursor(tbl, flds) as cursor:
    for lon, lat, x_wm, y_wm, x_nad, y_nad in cursor: 
        x_wm, y_wm = convertXY(lon,lat,latlon,webmer)
        x_nad, y_nad = convertXY(lon,lat,latlon,nad)
        cursor.updateRow((lon, lat, x_wm, y_wm, x_nad, y_nad))
del cursor

# using an insert cursor -and/or- print
coords = [ # can come from text file, etc.
    [175.5,-40.5],
    [135.5,-30.5],
    [90.5,-20.5],
    [45.5,-10.5],
    [0.0,0.0],
    [-45.5,10.5],
    [-90.5,20.5],
    [-135.5,30.5],
    [-175.5,40.5],
    ]
print "LON\tLAT\tX_WM\tY_WM\tX_NAD\tY_NAD"
cursor = arcpy.da.InsertCursor(tbl, flds)
for lon, lat in coords:
    x_wm, y_wm = convertXY(lon,lat,latlon,webmer)
    x_nad, y_nad = convertXY(lon,lat,latlon,nad)
    print "{}\t{}\t{}\t{}\t{}\t{}".format(lon, lat, x_wm, y_wm, x_nad, y_nad)
    cursor.insertRow((lon, lat, x_wm, y_wm, x_nad, y_nad))
del cursor‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
JoeBorgione
MVP Esteemed Contributor

Most excellant Randy Burton‌ !

That should just about do it....
0 Kudos
JoshuaBixby
MVP Esteemed Contributor

Randy Burton‌, looking at the code for convertXY, there are some risks with the structure of error handling.  Setting aside the larger debate on whether one should trap errors and hide, or partially hide, the specifics from users, there is a fairly easy way to give the user more information about what went wrong while still allowing the operation to continue.

With the function as written, there are at least 3 things that can generate errors (creating the Point, creating the PointGeometry, and projecting the PointGeometry), but the user will have no idea which step had the issue.  I suggest using information from the error itself in what is printed:

def convertXY(x, y, inSR, outSR):
    try:
        ptGeo = arcpy.PointGeometry(arcpy.Point(x,y),inSR).projectAs(outSR)
        return (ptGeo.centroid.X, ptGeo.centroid.Y)
    except BaseException as e:
        print(e)
        return (0.0,0.0)

On a much smaller note, since you are asking the user for x,y coordinates and not an ArcPy Point, what about asking for the EPSG codes instead of an ArcPy SpatialReference object. 

def convertXY(x, y, in_epsg, out_epsg):
    try:
        inSR = arcpy.SpatialReference(in_epsg)
        outSR =  arcpy.SpatialReference(out_epsg)
        ptGeo = arcpy.PointGeometry(arcpy.Point(x,y),inSR).projectAs(outSR)
        return (ptGeo.centroid.X, ptGeo.centroid.Y)
    except BaseException as e:
        print(e)
        return (0.0,0.0)
RandyBurton
MVP Regular Contributor


Great suggestions.  Thank you.  So far, the typical cause for an error in my testing was swapping the x and y coordinates.

0 Kudos