Select to view content in your preferred language

Convert lat-long to WKT spatial reference?

12535
20
Jump to solution
02-15-2016 02:34 PM
RickThiel
Regular 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
1 Solution

Accepted Solutions
RandyBurton
MVP Alum

Does something like this help:

import arcpy
# WGS 1984 : (4326) Lat/Lon
x = -89.384171
y = 43.074743
# xy = Wisconsin Capitol

sr = arcpy.SpatialReference(r'C:\path2file\custom.prj') # your projection file

ptGeometry = arcpy.PointGeometry(arcpy.Point(x,y),arcpy.SpatialReference(4326)).projectAs(sr)

print ptGeometry.JSON
print x, y
print ptGeometry.firstPoint.X, ptGeometry.firstPoint.Y
# converted xy = 2538768.50758 1604492.87376

View solution in original post

20 Replies
DanPatterson_Retired
MVP Emeritus

you still have to project from your initial coordinate system to your desired on...

Project—Help | ArcGIS for Desktop  

see the code snippet

However in the geometry classes, WKT etc are read-only properties

Geometry—Help | ArcGIS for Desktop 

but once projected, you can get it back out to JSON

but a simpler solution would be just to use the PointGeometry class and snag the point coordinates after  they are projected if all you need to know is the values

RickThiel
Regular Contributor

Thanks Dan for pointing me in the right direction. The example on the ESRI help page shows that when they Project a point, they immediately save it to an output feature class.  Is there any way to project the point feature on the fly without having to save it to an output file/db?

For example, the geo-coder returns a point with coordinates x=-90.079346, y=42.718911  (spatial ref = 4326).  Could I take those coordinates and project them on the fly to return the true coordinates that I want?

0 Kudos
DanPatterson_Retired
MVP Emeritus

of course

in the PointGeometry class

PointGeometry—Help | ArcGIS for Desktop

a PointGeometry object needs a point (long story), so you are there.

projectAs (spatial_reference, {transformation_name})

Projects a geometry and optionally applies a geotransformation.

To project, the geometry needs to have a spatial reference, and not have an UnknownCoordinateSystem. The new spatial reference system passed to the method defines the output coordinate system. If either spatial reference is unknown the coordinates will not be changed. The Z- and measure values are not changed by the ProjectAs method.

0 Kudos
RickThiel
Regular Contributor

OK, I think I understand.  But, because I am using a custom projection, won't that be a problem?  In the projectAs() method, it says "To project, the geometry needs to have a spatial reference, and not have an UnknownCoordinateSystem"

0 Kudos
DanPatterson_Retired
MVP Emeritus

but you do have a spatial reference in your example from the geocoder...  4326 by the coordinates.  That I presume is what they are coming in as...define it as such, then set the output projection

RebeccaStrauch__GISP
MVP Emeritus

Since I've been working on something similar the past couple days (not the json part), I'll throw this out there.  I'm not sure how you parsed your lat, longs, but this assuming a PtID, long, lat like

ptList = [[104, -148.631196, 64.010037], [102, -147.707779, 63.680584], [103, -147.469019, 63.975227]]

This script will take a structure like that, create a FC in DD, add the x/y as fields for kicks, then project to you desired wkid (mine is 3338, AK Albers).  I have also used AlterField (not shown here) on the field POINT_X, POINT_Y so I could use the addxy again in the projected FC (and a bunch of other fields).  If you don't have a PtID field,m or you input is a bit different, that will have to be modified.

See if this will help.

import os
import arcpy
from arcpy import env

theWorkspace = r'C:\__temp\test.gdb'
arcpy.env.workspace = theWorkspace   #not until final save
arcpy.env.overwriteOutput = True    

# pts with a ID, long, lat
ptList = [[104, -148.631196, 64.010037], [102, -147.707779, 63.680584], [103, -147.469019, 63.975227]]

ptFields = [["PtID", "SHORT", "5", "", ""]]
newFcs = [["newpts", ptFields]]    # , ["tmppt2", pt2Fields], ["tmppt3", pt2Fields]]
newFCMySR = "newPtsMySR"

mySR = arcpy.SpatialReference(3338)
geoSR = arcpy.SpatialReference(4269)

for newFC, fields in newFcs:
  print(newFC)
  if not arcpy.Exists(arcpy.os.path.join(theWorkspace, newFC)):
    print("{0} does not exists, creating....".format(newFC))
    arcpy.CreateFeatureclass_management(theWorkspace, newFC, "POINT", "", "", "", initialSR) 
    for addFld in fields:
      arcpy.AddField_management(newFC, addFld[0], addFld[1], addFld[2], addFld[3], addFld[4])
      print(" ..added field {0} to {1}".format(addFld[0], newFC))
    #listNewFields = arcpy.ListFields(newFC)
  else:
    print("{0} exists".format(newFC))

for newFC, newFields in newFcs:
  fldNames = []
  iCurFields = ['SHAPE@XY']
  for fn in newFields:
    print(fn[0])
    fldNames.append(fn[0])
    iCurFields.append(fn[0])
  print("\nField names in lists {0}: {1}".format(newFC,fldNames))
  print("   and for the cursor: {0}".format(iCurFields))
  cursor = arcpy.da.InsertCursor(newFC, iCurFields)
  print(" ->working with {0}".format(newFC))

  for ptID, x, y in ptList:
    print(ptID)
    xy = (x, y)
    print(xy)
    cursor.insertRow([xy, ptID])
  
  arcpy.AddXY_management(newFC)

arcpy.Project_management(newFC, newFCMySR, mySR)
RickThiel
Regular Contributor

Thanks Rebecca,

I started going down the same route you chose, by creating a temp GDB to store the new projected coordinates.  Late yesterday, I finally got a good result!  Yeah!  I was going to reply to your suggestion, but I was too tired last night and I went home instead.

Now, I see the suggestion from your fellow Alaskan, Randy Burton​ this morning and it seems like a much cleaner way to go. I think that is what I will use.  Time to rip out some code... 

Thanks Dan PattersonRebecca Strauch, GISP​ and Randy Burton​ for all of your help and attention.

Take care, --Rick

0 Kudos
RandyBurton
MVP Alum

Does something like this help:

import arcpy
# WGS 1984 : (4326) Lat/Lon
x = -89.384171
y = 43.074743
# xy = Wisconsin Capitol

sr = arcpy.SpatialReference(r'C:\path2file\custom.prj') # your projection file

ptGeometry = arcpy.PointGeometry(arcpy.Point(x,y),arcpy.SpatialReference(4326)).projectAs(sr)

print ptGeometry.JSON
print x, y
print ptGeometry.firstPoint.X, ptGeometry.firstPoint.Y
# converted xy = 2538768.50758 1604492.87376
RebeccaStrauch__GISP
MVP Emeritus

Randy,  thanks for that.  I think I'll go back to my process and see if I can use that instead of the projection process I've been doing.  Only reason I was doing that was to get the coords.  You way is nice and clean.

0 Kudos