Convert lat-long to WKT spatial reference?

11496
20
Jump to solution
02-15-2016 02:34 PM
RickThiel
Occasional Contributor

Hello,  arcpy newbie here... sorry if I am asking a basic question.  I am looking for help on transforming lat/long decimal degrees to a different spatial reference, specifically a well-known-text (WKT) custom reference that we are using.

I have sent a street address to an address locator (findAddressCandidates) and it returns JSON back to me.  If the candidate score is near 100, then I parse the JSON to grab the lat/long.  That all works fine.  But, now I would like to use an insertCursor to add this point to my feature class.  I am unclear how I would convert a lat long to the custom spatial ref., WKT?

I have been trying many different things in the past few days but nothing works for me.  I have tried passing the WKT in the URL for the address locator.  It seems like the locator does not recognize the WKT and converts it to the default 4326 spatial reference (not what I want).

FYI.. this is the custom spatial reference that I am using:

PROJCS["NAD_1983_Lambert_Conformal_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",1840330.965],PARAMETER["False_Northing",527997.6316666666],PARAMETER["Central_Meridian",-92.0],PARAMETER["Standard_Parallel_1",41.15],PARAMETER["Standard_Parallel_2",45.0],PARAMETER["Latitude_Of_Origin",40.15],UNIT["Foot_US",0.3048006096012192]]

Thanks, --Rick

0 Kudos
20 Replies
RandyBurton
MVP Alum

If you know the reference numbers, you can use this instead:

# WGS 1984 : (4326) Lat/Lon
# WGS 1984 Web Mercator (auxiliary sphere) : (102100) or (3857)
ptGeometry = arcpy.PointGeometry(arcpy.Point(x,y),arcpy.SpatialReference(4326)).projectAs(arcpy.SpatialReference(3857))
0 Kudos
RickThiel
Occasional Contributor

Thanks Randy!  That's exactly what I was looking for!

Nice touch using the Wisconsin Capitol.    Rebecca's lat/long examples were SO confusing to me!  Haha.  I mean, -147.707779, 63.680584, what is that???  Sara Palin's house?

Take care everybody!  Thanks again!  -_Rick

0 Kudos
JohnMcoy
New Contributor III

hi Randy, how can I take x,y coordinates from featureclass? Because I don`t have xy values in code.

0 Kudos
RandyBurton
MVP Alum

Using a search cursor, you can loop through the points in a feature class, read the geometry and convert to longitude/latitude.  Something like this:

for row in arcpy.da.SearchCursor("MyPointLayer",["OID@","SHAPE@XY"]):

     print("Feature {}:".format(row[0])),

     px, py = row[1]
     print("{}, {}".format(px, py)),

     # WGS 1984 : (4326) Lat/Lon  
     # WGS 1984 Web Mercator (auxiliary sphere) : (102100) or (3857)  
     ptGeometry = arcpy.PointGeometry(arcpy.Point(px,py),arcpy.SpatialReference(3857)).projectAs(arcpy.SpatialReference(4326))

     print ptGeometry.firstPoint.X, ptGeometry.firstPoint.Y‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

This code may be helpful in getting the spatial reference information for the point feature:

SR = arcpy.Describe("MyPointLayer").spatialReference
print SR.name
print SR.factoryCode‍‍‍
JohnMcoy
New Contributor III

Could you explain, how to separate coordinates in different columns ? Now I have coordinates like this(photo). I want 4 different columns Point X Point Y in WGS1984 and POINT X POINT Y in WGS 1984 Web Mercator..

0 Kudos
RandyBurton
MVP Alum

Here's a version of the previous script that puts the x and y coordinates into a comma delimited format.

# print a header row - comma delimited
print "{}, {}, {}, {}, {}".format(
    '"OBJECTID"',
    '"X_WebM"',
    '"Y_WebM"',
    '"X_Lon"',
    '"Y_Lat"'
    )

# we will assume the feature is in web mercator
for row in arcpy.da.SearchCursor("MyPointLayer",["OID@","SHAPE@XY"]):

    # SHAPE@XY contains both x and y coordinates
    px, py = row[1]

    # convert px and py to longitude and latitude
    # WGS 1984 : (4326) Lat/Lon  
    # WGS 1984 Web Mercator (auxiliary sphere) : (102100) or (3857)  
    ptGeometry = arcpy.PointGeometry(arcpy.Point(px,py),
                                     arcpy.SpatialReference(3857)).projectAs(arcpy.SpatialReference(4326))

    print "{}, {}, {}, {}, {}".format(
        row[0], # OBJECTID
        px, # X_WebM
        py, # Y_WebM
        ptGeometry.firstPoint.X, # X_Lon
        ptGeometry.firstPoint.Y # Y_Lat
        )

This should make it easier to see how to separate the geometry into columns.  You can add other attribute fields to the list in the Search Cursor:

for row in arcpy.da.SearchCursor("MyPointLayer",["OID@","SHAPE@XY","FieldName"]):
    ....
    row[2] # FieldName‍‍‍
DanPatterson_Retired
MVP Emeritus

Rick ...

if you want to weigh in here anytime... the suggestions are piling up... I will whip in a numpy solution if this isn't closed soon

0 Kudos
RebeccaStrauch__GISP
MVP Emeritus

Randy's solution was great for me (I do use the SR #'s) except since they were just one points I used .centroid instead of .startpoint.   this saved two projections/addXY/AlterField on 3 different files for me.    This is why I really like geonet.

If anyone has a simple solution using the geometry for converting the long, lat values to DDM and DMS (into new fields) other than using the Convert Coordinate Notation—Help | ArcGIS for Desktop    which works great and adds the new fields, but to a new FC each time (3 times x 3 files in my case), I would love to see it.  I may just do the conversion/formatting in a and add the values to my [ ]. 

Dan....simple Numpy suggestion perhaps? DD  (negative long in my case) -> DDM and DMS?  I can fire this off to a new question if needed, but I know there are already thread on the conversion (and I have snippets already)

0 Kudos
RandyBurton
MVP Alum

Getting off topic a little... Rebecca, you could write a function.  For DD to DMS, something like:

def dd2dms(dd,xy):
    is_positive = dd >= 0
    dd = abs(dd)
    m,s = divmod(dd*3600,60)
    d,m = divmod(m,60)
    if is_positive:
        if xy == 'lat' : nsew = ' N'
        else : nsew = ' E'
    else:
        if xy == 'lat' : nsew = ' S'
        else : nsew = ' W'
    return str(int(d))+' '+str(int(m))+' '+str(round(s,1))+ nsew

print 61.217176, -149.889006
print dd2dms(61.217176,'lat'), dd2dms(-149.889006,'lon')
# 61 13 1.8 N 149 53 20.4 W
0 Kudos
DanPatterson_Retired
MVP Emeritus

Rebecca, you probably haven't had time to scroll around the Numpy Repository, but I have this in there

Converting Decimal Degrees to DMS format (there is a zip as well)

0 Kudos