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.
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
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'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.
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.
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.
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
Most excellant Randy Burton !
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)
Great suggestions. Thank you. So far, the typical cause for an error in my testing was swapping the x and y coordinates.