topic Re: How do you project a single pointGeometry object with arcpy? in Python Questions

How do you project a single pointGeometry object with arcpy?

JohnMcoy — Sat, 11 Dec 2021 09:37:13 GMT

Hi everyone,

So I got a task, and I don`t know how to do it. In csv file, i have to create two new columns which shows WGS-84 coordinates. Maybe someone, have solution for this ? I can`t undertstand what is in 29 line... I tried to google it.. And I found the example... but without explanation.

Thanks for any help

import arcpy 
import csv

fc = r"#"
out_csv = r"#"

fields = {
    "OBJECTID": "OBJECTID",
    "Global_ID": "Global_ID",
    "Global_Copy": "Global_Copy",
    "SHAPE@X": "POINT_X",
    "SHAPE@Y": "POINT_Y",
    "SHAPE@X": "LAT",
    "SHAPE@Y": "LONG"
}

table_fields = []
csv_fields = []
for field in fields:
    table_fields.append(field)
    csv_fields.append(fields[field])

def tableToCSV(fc, out_csv):
    with open(out_csv, 'wb') as csv_file:
        writer = csv.writer(csv_file, delimiter=';', lineterminator='\n')
        writer.writerow(csv_fields)
        with arcpy.da.SearchCursor(fc, table_fields) as cursor:
            for row in cursor:
                arcpy.PointGeometry(arcpy.Point(X,Y),arcpy.SpatialReference(3346)).projectAs(arcpy.SpatialReference(4326))
                point = Point(row[idxLAT],row[idxLONG])
                #TODO: reproject
                reppoint = point
                row[idxLAT] = reppoint.x
                row[idxLONG] = reppoint.y
                
                writer.writerow(row)
        print out_csv  
    #csv_file.close()

tableToCSV(fc, out_csv)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: How do you project a single pointGeometry object with arcpy?

NeilAyres — Tue, 02 May 2017 07:24:45 GMT

Well WKID 3346 appears to be some coordinate system in Lithuania.

LKS94 / Lithuania TM: EPSG Projection -- Spatial Reference

Line 29, is taking 2 values (X, Y), creating a point geometry in 3346 and then re-projecting them to WKID 4326, which is GCS_WGS84.

If you already had a geometry object, the general form would be :

geom2 = geom1.projectAs(newSR)‍

Not sure in the original script where X & Y are coming from.

Re: How do you project a single pointGeometry object with arcpy?

NeilAyres — Tue, 02 May 2017 07:27:12 GMT

You could just import the csv, make a point feature out of it, then re-project it to GCS_WGS84 directly.

If you want a programmatic solution, then some details of what the csv file looks like would be good.

Re: How do you project a single pointGeometry object with arcpy?

JohnMcoy — Tue, 02 May 2017 07:39:58 GMT

I have database with ~50 columns. I take about 20 specific columns. As you see in array-fields. In csv Poin-X and Point-Y has show 3346 cordinates and Lat, Long has show 4326. In CSV file must create two new columns Lat and Long. CSV looks like this i hope you understand everything. Is there enough information ??

Re: How do you project a single pointGeometry object with arcpy?

NeilAyres — Tue, 02 May 2017 07:49:03 GMT

So that csv is already exported from a feature class? It has an OBJECTID column.

And the presence of POINT_X, POINT_Y would also suggest that this has been used...

Add Geometry Attributes—Help | ArcGIS Desktop

So, if you haven't got the original feature instead of this exported version, why not just re-import it into ArcMap?

There are several ways of adding the XY (in whatever form) back to the point feature.

Generally, if the required coordinates are not the actual coordinates of the point feature, I use ArcMap, add the feature, change the coord sys of the data frame to the system I want, then use the calculate geometry option on 2 new columns, use the coord sys of the data frame.

Re: How do you project a single pointGeometry object with arcpy?

JohnMcoy — Tue, 02 May 2017 07:59:06 GMT

No, the script has export csv file with all needed columns. X and Y coordinates coming from featureclass with SHAPE@X and SHAPE@Y. And I can`t understand, how to take the points coordinates from fc and put in this place... http://prntscr.com/f35bv3

Re: How do you project a single pointGeometry object with arcpy?

NeilAyres — Sat, 11 Dec 2021 09:37:16 GMT

Okay, so lets have a look at that script so far...

The definition of the "fields" dictionary isn't making a lot of sense.

fields = {
    "OBJECTID": "OBJECTID",
    "Global_ID": "Global_ID",
    "Global_Copy": "Global_Copy",
    "SHAPE@X": "POINT_X",
    "SHAPE@Y": "POINT_Y",
    "SHAPE@X": "LAT",
    "SHAPE@Y": "LONG"
}‍‍‍‍‍‍‍‍‍

Each row is a pair of key : value. Supposedly to map fields in the Fc with output fields in the csv. But why would "SHAPE@X/Y" appear twice? Because rows come out of a dict in an arbitrary order, what is mapping the XY to the name in the csv?

The tokens "SHAPE@X" & "SHAPE@Y" recover the X & Y from the geometry column of the Fc and will be whatever coordinate system the Fc is in.

Reading geometries—Help | ArcGIS Desktop

So, maybe this little example will give you an idea. I made a small point feature (5 points), in UTM35S_WGS84.

I can read it and print the XY metres, and also the Long/Lat.

import sys, os, arcpy

inFgdb = r"C:\Data\ESRI-SA\ForumProblem\scratch.gdb"
inFc = "TestPoints_UTM"

inData = os.path.join(inFgdb, inFc)

inFlds = ["OID@", "SHAPE@"]

srIN = arcpy.Describe(inData).SpatialReference
print "Input SR is {}".format(srIN.name)

srOUT = arcpy.SpatialReference(4326) # GCS_WGS84
print "Output SR is {}".format(srOUT.name)


with arcpy.da.SearchCursor(inData, inFlds) as Cur:
    for row in Cur:
        ID = row[0] # first item of inFlds list
        geom = row[1] # second item of inFlds list
        # we have used the geometry object, not the XY tokens
        X, Y = geom.centroid.X, geom.centroid.Y
        geom2 = geom.projectAs(srOUT)
        X2, Y2 = geom2.centroid.X, geom2.centroid.Y
        # print out the coordinates, formatted
        # 1 decimal for the metres based UTM coords
        # 6 decimals for the degree based coords
        print "ID {}, UTM metres {:.1f} {:.1f}, Long/Lat {:.6f} {:.6f}".format(ID, X, Y, X2, Y2)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Gives me this output

>>> 
Input SR is WGS_1984_UTM_Zone_35S
Output SR is GCS_WGS_1984
ID 1, UTM metres 542579.7 7064795.0, Long/Lat 27.427425 -26.536807
ID 2, UTM metres 562671.5 7084560.1, Long/Lat 27.628134 -26.357606
ID 3, UTM metres 590440.6 7078026.2, Long/Lat 27.906892 -26.415105
ID 4, UTM metres 595994.4 7060057.9, Long/Lat 27.963929 -26.576958
ID 5, UTM metres 574759.2 7034085.6, Long/Lat 27.752251 -26.812725
>>> ‍‍‍‍‍‍‍‍‍

In your script you are using the XY tokens, so have to construct a point geometry from them so it can be projected.

I chose do do it using the geometry object itself.

Hope this helps.

Re: How do you project a single pointGeometry object with arcpy?

JohnMcoy — Tue, 02 May 2017 14:10:15 GMT

Thanks for example Neil, but could you show me, how to separate ALL coordinates in different columns ? Thanks!

Re: How do you project a single pointGeometry object with arcpy?

NeilAyres — Tue, 02 May 2017 15:20:26 GMT

Sorry, but I am not quite with you here....

You now have a bunch of variables, you can construct a list of those values, then push it through the csv writer to write each line.

Re: How do you project a single pointGeometry object with arcpy?

RandyBurton — Sat, 11 Dec 2021 09:37:20 GMT

Hello John,

When I export a point feature to a CSV or Excel table, I like to have the X and Y coordinates in longitude and latitude. To do this, I use the Add Geometry Attributes tool before exporting my table. This tool adds two columns, Point_X and Point_Y, to your feature. By specifying a spatial reference, you can have the XY in longitude/latitude. This tool does modify the input table by adding these columns. You could make a copy of your feature, perhaps "in-memory", and modify it if you did not want the columns in your original feature.

I would use the following code in your script before your lines to export to csv to create and populate the feature with Point_X and Point_Y fields:

spatial_reference = arcpy.SpatialReference("WGS 1984")
arcpy.AddGeometryAttributes_management(fc,"POINT_X_Y_Z_M","#","#",spatial_reference)‍‍‍‍‍‍‍‍

Like Neil Ayres‌, I am confused by your files dictionary. I am not sure why you are adding a second SHAPE@X/Y (X is Long and Y is Lat), so I would suggest dropping this.

Instead of using Add Geometry Attributes, you might try something like this to get longitude and latitude:

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‍‍‍‍‍‍‍‍‍‍‍‍

Hope this helps.

Re: How do you project a single pointGeometry object with arcpy?

RandyBurton — Sat, 11 Dec 2021 09:37:22 GMT

This provides some interesting results. By specifying a spatial reference in the search cursor, you can get latitude and longitude very easily. And it is easy to have tab delimited output.

# WGS 1984 : (4326) Lat/Lon  

print "{}\t{}\t{}".format(
    "OBJECTID",
    "POINT_X",
    "POINT_Y"
    )
for row in arcpy.da.SearchCursor("MyPointLayer",["OID@","SHAPE@X","SHAPE@Y"],spatial_reference=arcpy.SpatialReference(4326)):

    print "{}\t{}\t{}".format(
        row[0], # OBJECTID
        row[1], # POINT_X (LON)
        row[2]  # POINT_Y (LAT)
    )
‍‍‍‍‍‍‍‍‍‍‍‍‍‍