topic Re: projectAs help in Python Questions

projectAs help

JamesCrandall — Mon, 13 May 2013 12:56:06 GMT

I'm attempting to implement some of the functionality found here:

http://support.esri.com/en/knowledgebase/techarticles/detail/40730

I need to acquire the x/y values where mouse click occurs, convert the values into Latitude/Longitude and then report these values.

The plan is to require a layer to be loaded in the map TOC, so there will be a defined spatial reference to go from. No problem getting these projected x/y values so far but I am not able to project these to WGS1984. Is this even possible? Or will I have to use in_memory in some capacity here?

      def onMouseDownMap(self, x, y, button, shift):         print "x: " + str(x) + " y: " + str(y) ## this reports as expected!                          ## get a layer loaded in the TOC to determine the spatial reference                 mxd = arcpy.mapping.MapDocument("CURRENT")         df = arcpy.mapping.ListDataFrames(mxd)[0]         testCount = len(arcpy.mapping.ListLayers(mxd, "", df))         if testCount==0:             msg = "There are no features selected to report"             pythonaddins.MessageBox(msg, 'Report Latitude/Longitude Values', 0)         else:               for lyr in arcpy.mapping.ListLayers(mxd):             dsc = arcpy.Describe(lyr)             spref = dsc.spatialReference              lpoint = arcpy.Point()             lpoint.X = x             lpoint.Y = y             ptGeometry = arcpy.PointGeometry(lpoint)             #ptGeom = arcpy.PointGeometery(arcpy.point(x,y),spref, False, False)              to_sr = arcpy.SpatialReference('WGS 1984')                                  projectedPoint = ptGeometry.projectAs(to_sr, r'NAD_1983_HARN_To_WGS_1984')             prjX = projectedPoint.X             prjY = projectedPoint.Y             print "X:" + str(prjX) + " Y:" + str(prjY)

The code above fails on the line:

[INDENT]projectedPoint = ptGeometry.projectAs(to_sr, r'NAD_1983_HARN_To_WGS_1984')[/INDENT]

Error msg:
"projectedPoint = ptGeometry.projectAs(to_sr, r'NAD_1983_HARN_To_WGS_1984')
File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 761, in projectAs
return convertArcObjectToPythonObject(self._arc_object.ProjectAs(*gp_fixargs((spatial_reference, transformation_name))))
ValueError: NAD_1983_HARN_To_WGS_1984"

Re: projectAs help

TimDine — Sat, 11 Dec 2021 14:09:14 GMT

Check that you're using a valid transformation for the projections you're using by feeding the shapefile to the project tool in a model and exporting it. I tried two variations of a script with different options for the transformation parameter. The transformation you're using blew up with the same error you're getting, and a transformation that I used from the exported script worked perfectly. I got the error you're using when I had a typo in my transformation string as well.

This code blows up with your error

import arcpy
theFC = r"c:\temp\myFeatureclass.shp"
g = arcpy.Geometry()
to_sr = arcpy.SpatialReference('WGS 1984')
geometryList = arcpy.CopyFeatures_management(theFC, g)
for geometry in geometryList:
 projectedGeometry = geometry.projectAs(to_sr,r'NAD_1983_HARN_To_WGS_1984')
 print projectedGeometry.centroid.X, projectedGeometry.centroid.Y

This code works great

import arcpy
theFC = r"c:\temp\myFeatureclass.shp"
g = arcpy.Geometry()
to_sr = arcpy.SpatialReference('WGS 1984')
geometryList = arcpy.CopyFeatures_management(theFC, g)
for geometry in geometryList:
 projectedGeometry = geometry.projectAs(to_sr,'NAD_1983_To_NAD_1983_CSRS_4 + NAD_1983_CSRS_To_WGS_1984_2')
 print projectedGeometry.centroid.X, projectedGeometry.centroid.Y

Re: projectAs help

JamesCrandall — Sat, 11 Dec 2021 14:09:17 GMT

Thanks, tim. However you are still requiring an input shapefile, but I will go back and see if I can rework things prior to my new approach.

I was hoping to avoid this and simply convert the geometry to a different spatial reference (WGS1984). Doesn't look like this is possible, so I came up with an "in_memory" implementation that seems to do what I need. The first click event takes about 3 seconds to build the in_memory FeatureClass but subsequent clicks the Lat/Lon values are reported immediately.

Solution:


    def onMouseDownMap(self, x, y, button, shift):
        print "x: " + str(x) + " y: " + str(y)
        
        ## get a layer loaded in the TOC to determine the spatial reference        
        mxd = arcpy.mapping.MapDocument("CURRENT")
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        testCount = len(arcpy.mapping.ListLayers(mxd, "", df))
        if testCount==0:
            msg = "There are no features selected to report"
            pythonaddins.MessageBox(msg, 'Report Latitude/Longitude Values', 0)
        else:    
          for lyr in arcpy.mapping.ListLayers(mxd):
            dsc = arcpy.Describe(lyr)
            spref = dsc.spatialReference
            
            lpoint = arcpy.Point()
            lpoint.X = x
            lpoint.Y = y
            ptGeometry = arcpy.PointGeometry(lpoint, spref)

            ptFC = "in_memory\lpoint"
            if arcpy.Exists(ptFC):
                arcpy.Delete_management(ptFC)
                
            arcpy.CreateFeatureclass_management("in_memory", "lpoint", "POINT", '', "DISABLED", "DISABLED",  spref, '')
            
            cursor = arcpy.da.InsertCursor(ptFC, "SHAPE@XY")
            cursor.insertRow((lpoint,))
            del cursor
            
            to_sr = r'GEOGCS["GCS_WGS_1984",' + \
                     'DATUM["D_WGS_1984",' + \
                     'SPHEROID["WGS_1984",6378137,298.257223563]],' + \
                     'PRIMEM["Greenwich",0],' + \
                     'UNIT["Degree",0.017453292519943295]]'
           
            flds = ["SHAPE@X", "SHAPE@Y"]
            scursor = arcpy.da.SearchCursor(ptFC, flds, '', to_sr)
            prjx = 0.00
            prjy = 0.00
            for row in scursor:
                #print "X: " + str(row[0]) + " Y: " + str(row[1])
                prjx = row[0]
                prjy = row[1]

            del scursor

            
            msg = "Latitude: " + str(prjx) + "\n" + "Longitude: " + str(prjy)
            pythonaddins.MessageBox(msg, 'Report Latitude/Longitude Values', 0)

this appears to be a good solution, but if anyone can spot any issue with the above, I'd appreciate it.

Re: projectAs help

TimDine — Mon, 13 May 2013 16:11:13 GMT

The fact an input shapefile is being used shouldn't be relevant. As long as you get a pointgeometry object somehow with a defined spatial reference the projectAs should work fine.

Your first solution shouldn't have any issue other than something being wrong with your transformation parameter (r'NAD_1983_HARN_To_WGS_1984'). I get your same error if I use the shapefile and your transformation. The samples I provided were simplified to use shapefiles to eleminate the rest of the code as being an issue.

Re: projectAs help

JamesCrandall — Mon, 13 May 2013 16:40:35 GMT

The fact an input shapefile is being used shouldn't be relevant. As long as you get a pointgeometry object somehow with a defined spatial reference the projectAs should work fine.

Your first solution shouldn't have any issue other than something being wrong with your transformation parameter (r'NAD_1983_HARN_To_WGS_1984'). I get your same error if I use the shapefile and your transformation. The samples I provided were simplified to use shapefiles to eleminate the rest of the code as being an issue.

Understood -- I will re-work the original attempt.

Thanks again,
j

Re: projectAs help

ChrisSnyder — Mon, 13 May 2013 16:46:23 GMT

It makes the code a bit harder to comprehend, but how about using the WKID (well known id) for your projections and transformation definitions? I have notices some "discrepencies" in the naming conventions, and well, computers like numbers, right? I just pull the codes from these documents:

C:\Program Files (x86)\ArcGIS\Desktop10.1\Documentation\geographic_transformations.pdf
C:\Program Files (x86)\ArcGIS\Desktop10.1\Documentation\geographic_coordinate_systems.pdf
C:\Program Files (x86)\ArcGIS\Desktop10.1\Documentation\projected_coordinate_systems.pdf

Re: projectAs help

JamesCrandall — Sat, 11 Dec 2021 14:09:20 GMT

The fact an input shapefile is being used shouldn't be relevant. As long as you get a pointgeometry object somehow with a defined spatial reference the projectAs should work fine.

Your first solution shouldn't have any issue other than something being wrong with your transformation parameter (r'NAD_1983_HARN_To_WGS_1984'). I get your same error if I use the shapefile and your transformation. The samples I provided were simplified to use shapefiles to eleminate the rest of the code as being an issue.

Tim -- I tried the code below but getting the same error for the transformation:

line 217, in onMouseDownMap
projectedPoint = ptGeometry.projectAs(to_sr, 'NAD_1983_To_NAD_1983_CSRS_4 + NAD_1983_CSRS_To_WGS_1984_2')
File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 761, in projectAs
return convertArcObjectToPythonObject(self._arc_object.ProjectAs(*gp_fixargs((spatial_reference, transformation_name))))
ValueError: NAD_1983_To_NAD_1983_CSRS_4 + NAD_1983_CSRS_To_WGS_1984_2

Code:


   lyr = arcpy.mapping.ListLayers(mxd)[0]
   dsc = arcpy.Describe(lyr)
   spref = dsc.spatialReference
            
   lpoint = arcpy.Point()
   lpoint.X = x
   lpoint.Y = y
   ptGeometry = arcpy.PointGeometry(lpoint, spref)
   to_sr = arcpy.SpatialReference('WGS 1984')
     
   ## fails on this line               
   projectedPoint = ptGeometry.projectAs(to_sr, 'NAD_1983_To_NAD_1983_CSRS_4 + NAD_1983_CSRS_To_WGS_1984_2')
   
   prjX = projectedPoint.X
   prjY = projectedPoint.Y
   print "X:" + str(prjX) + " Y:" + str(prjY)

Re: projectAs help

JamesCrandall — Sat, 11 Dec 2021 14:09:22 GMT

It makes the code a bit harder to comprehend, but how about using the WKID (well known id) for your projections and transformation definitions? I have notices some "discrepencies" in the naming conventions, and well, computers like numbers, right? I just pull the codes from these documents:

C:\Program Files (x86)\ArcGIS\Desktop10.1\Documentation\geographic_transformations.pdf
C:\Program Files (x86)\ArcGIS\Desktop10.1\Documentation\geographic_coordinate_systems.pdf
C:\Program Files (x86)\ArcGIS\Desktop10.1\Documentation\projected_coordinate_systems.pdf

Thanks for your input!

I attempted to set a WKID value but end up with the same error:

projectedPoint = ptGeometry.projectAs(to_sr, r'1900')
File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 761, in projectAs
return convertArcObjectToPythonObject(self._arc_object.ProjectAs(*gp_fixargs((spatial_reference, transformation_name))))
ValueError: 1900

            lyr = arcpy.mapping.ListLayers(mxd)[0]
            dsc = arcpy.Describe(lyr)
            spref = dsc.spatialReference
            
            lpoint = arcpy.Point()
            lpoint.X = x
            lpoint.Y = y
            ptGeometry = arcpy.PointGeometry(lpoint, spref)
            to_sr = arcpy.SpatialReference('WGS 1984')
                    
            projectedPoint = ptGeometry.projectAs(to_sr, 1900) #<--attempted with both integer value and r'1900'
            prjX = projectedPoint.X
            prjY = projectedPoint.Y
            print "X:" + str(prjX) + " Y:" + str(prjY)

The from system is NAD83 HARN Florida East FIPS 901 Feet, so perhaps my ignorance on the details for these coordinate conversions is just completely the cause here.

Re: projectAs help

ChrisSnyder — Mon, 13 May 2013 17:09:38 GMT

Just use the integer (don't format as string):

projectedPoint = ptGeometry.projectAs(12345, 54321)

Re: projectAs help

JamesCrandall — Sat, 11 Dec 2021 14:09:25 GMT

Just use the integer (don't format as string):
projectedPoint = ptGeometry.projectAs(12345, 54321)

Thanks again...

From the geographic_coordinate_systems.pdf (pg 35) the WKID for "WGS_1984" = 7030
From the geographic_transformations.pdf (pg 66) the WKID for "NAD_1983_HARN_To_WGS_1984_2" = 1900

Here's the updated code only using the WKID's:


            lyr = arcpy.mapping.ListLayers(mxd)[0]
            dsc = arcpy.Describe(lyr)
            spref = dsc.spatialReference
            
            lpoint = arcpy.Point()
            lpoint.X = x
            lpoint.Y = y
            ptGeometry = arcpy.PointGeometry(lpoint, spref)
                   
            projectedPoint = ptGeometry.projectAs(7030, 1900)
            prjX = projectedPoint.X
            prjY = projectedPoint.Y
            print "X:" + str(prjX) + " Y:" + str(prjY)

Here's the output error:

projectedPoint = ptGeometry.projectAs(7030, 1900)
File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 761, in projectAs
return convertArcObjectToPythonObject(self._arc_object.ProjectAs(*gp_fixargs((spatial_reference, transformation_name))))
RuntimeError: Object: CreateObject error creating spatial reference

Re: projectAs help

ChrisSnyder — Mon, 13 May 2013 17:53:43 GMT

I think you want GCS_WGS_1984 = 4326 (on page 18). I messed around with it and it seems to work for me like this:

pntObj = arcpy.Point(X= 12345,Y=6789)  inSr = arcpy.SpatialReference(2927) #NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet outSr = arcpy.SpatialReference(4326) #GCS_WGS_1984 pntGeom = arcpy.PointGeometry(pntObj, inSr) wgs84PntGeom = pntGeom.projectAs(outSr, "NAD_1983_HARN_To_WGS_1984_2")

Which yields:

>>> wgs84PntGeom.JSON
u'{"x":-126.81813166550971,"y":45.173061523748544,"spatialReference":{"wkid":4326}}'

Re: projectAs help

JamesCrandall — Sat, 11 Dec 2021 14:09:28 GMT

I think you want GCS_WGS_1984 = 4326 (on page 18). I messed around with it and it seems to work for me like this:
pntObj = arcpy.Point(X= 12345,Y=6789) 
inSr = arcpy.SpatialReference(2927) #NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet
outSr = arcpy.SpatialReference(4326) #GCS_WGS_1984
pntGeom = arcpy.PointGeometry(pntObj, inSr)
wgs84PntGeom = pntGeom.projectAs(outSr, "NAD_1983_HARN_To_WGS_1984_2")
Which yields:

>>> wgs84PntGeom.JSON
u'{"x":-126.81813166550971,"y":45.173061523748544,"spatialReference":{"wkid":4326}}'

Thanks for this -- once I figured out my inSr WKID, it returns the lat/lon correctly. I went ahead and made your post the answer, but...

Can you suggest how to parse out the actual x/y values out of that JSON string?

Thanks again for your help.
j

Re: projectAs help

JamesCrandall — Sat, 11 Dec 2021 14:09:30 GMT

Can you suggest how to parse out the actual x/y values out of that JSON string?

Pretty simple to do:


      data = json.loads(wgs84PntGeom.JSON)
      wgsX = data["x"]
      wgsY = data["y"]
      print "X:" + str(wgsX) + " Y:" + str(wgsY)

Re: projectAs help

ChrisSnyder — Sat, 11 Dec 2021 14:09:33 GMT

I think this would be the "best" way. Since it's a point you could also use the .centroid property (and some other ones too).

partObj = wgs84PntGeom.getPart(0)
partObj.X
partObj.Y

So, just to clarify - probably best not to parse the JSON... I just put that in there to illustrate it was returning WGS84 coordinates.

Re: projectAs help

JamesCrandall — Sat, 11 Dec 2021 14:09:36 GMT

I think this would be the "best" way. Since it's a point you could also use the .centroid property (and some other ones too).
partObj = wgs84PntGeom.getPart(0)
partObj.X
partObj.Y
So, just to clarify - probably best not to parse the JSON... I just put that in there to illustrate it was returning WGS84 coordinates.

Understood. Although that JSON string is fast enough to print out for the onMouseMoveMap event as the user moves across the view!

Thanks again