Calc Latitude, Longitude Field

[BlakeTerhune]

After getting a lot of requests to get latitude and longitude coordinates for features, I created a generic script that creates two new fields in the existing feature feature class and and calcs the latitude and longitude of the feature's "true centroid." The only problem is that you would have to run this script any time a feature's geometry was changed. I did some searching for something that already existed and didn't find anything that was simple and generic enough, so please let me know if you can do it better.

import arcpy

fc = r"C:\temp\working.gdb\temp_fc"

# Define and create longitude/latitude fields
## Longitude is X, Latitude is Y
fieldsNew = ["GCS_WGS_1984_X", "GCS_WGS_1984_Y"]
fieldsExisting = {f.name: f.type for f in arcpy.ListFields(fc)}
for fNew in fieldsNew:
    # Add field if it doens't already exist
    if fNew not in fieldsExisting.keys():
        arcpy.AddField_management(
            fc,  ## in_table
            fNew,  ## field_name
            "DOUBLE",  ## field_type
        )
    # If data type of existing field is not a Double, throw Exception
    elif fieldsExisting[fNew] != "Double":
        raise Exception(
            "Existing {} field is type {}. Must be 'Double'".format(
                fNew, fieldsExisting[fNew]
            )
        )
# Update fields with GCS_WGS_1984 lat/long
print("Updating Field values")
fieldsNew.append("SHAPE@")
GCS_WGS_1984 = arcpy.SpatialReference(4326)
with arcpy.da.UpdateCursor(fc, fieldsNew, "", GCS_WGS_1984) as u_cursor:
    for row in u_cursor:
        point = row[2].trueCentroid
        row[0] = point.X
        row[1] = point.Y
        u_cursor.updateRow(row)

python snippets​

[XanderBakker]

Did you look at: Add Geometry Attributes (Data Management)?

Question: what if your input featureclass does not have a sr WGS 1984 and it needs a specific transformation? Will the output coordinates be correct?

For polylines I prefer to use the midpoint of the line (positionAlongLine (0.5, True)) over the center of gravity...

[BlakeTerhune]

Great points, Xander. I don't think I saw AddGeometryAttributes_management(). Or maybe I did but didn't see that you could specify the coordinate system. I'll do some testing with that!

I have not needed to use a a specific transformation so have not tested it. What is a scenario I could use to test this? I don't understand what you mean when you say it "does not have a sr WGS 1984."

For the lines, I agree, a midpoint would be better. I just didn't want to write extra code to check the geometry type and do different things. Point features are obviously preferred so anything else requires some compromise. That would be a good enhancement though, to do different calculations of feature to point conversion based on the geometry type.

I also thought maybe instead of hard-coding the field name, I could just make a new field using the current date/time in the name. That would let me get rid of the field name and data type check and it would also hint that the field's values were calculated in a batch process separately from editing the geometry.

[NeilAyres]

If the datum (underlying GCS of the coordinate system) of the input is not WGS84, then a transformation is required if you want WGS84 based coordinates to come out.

I have never actually used the spatial reference option when opening an Update cursor. Generally, the feature exists and it knows what coordinate system it is in.... The help files here are unclear as to how this option affects the outcome.

I have in the past used the .projectAs(SR2) method on the geometry object to re-project geometries on the fly inside python. Specifying a transform is a little trickier and you have to set this via environment setting.

[XanderBakker]

Neil Ayres is right that you can project a geometry using you code and apply the transformation to it. However, it is not necessary to do this in the environment settings as you can use the the second (optional) parameter of the arcpy.Geometry::projectAs (spatial_reference, {transformation_name}) function.

An example of requiring a transformation is projecting data in Colombia from GSC_Bogota to MAGNA_Colombia_Bogota. This would require multiple steps (projections) and when going from GCS_MAGNA to GCS_Bogota, the transformation Bogota_To_MAGNA_Region_5_MB should be applied (according to specifications here).

If the transformation is not specified another transformation could be applied resulting in an inaccurate on-the-fly projection of the geometry and yielding inaccurate output coordinates.

In case these types of problems do not affect you, there is no problem. However, specifying a different spatial reference (sr) then your input as parameter of a cursor may have unexpected side-effects.

About the field names, you could use a very specific name that a normal user would unlikely create him/herself or create a tool and let the user define the names of the output fields. In case they exist verify it during validation or simply return an error.

[BlakeTerhune]

Is there an easier way to just convert state plane coordinates to decimal degrees without reprojecting? Our data is in NAD_1983_HARN_StatePlane_Arizona_Central_FIPS_0202.

EDIT:

Looks like the accepted transformation from NAD_1983 to WGS_1984 in the contiguous 48 states of the US is NAD_1983_To_WGS_1984_5. However, I don't see how to specify a transformation when using AddGeometryAttributes_management().

[XanderBakker]

If the tool itself does not provide it, you can use the environment settings:

• Output Coordinate System (Environment setting)
• Geographic Transformations (Environment setting)

[DanPatterson_Retired]

Does this tool not modify a file in-place?  That could be the issue.  in memory appears to be used for outputs unless the inputs are there as there as well

[BlakeTerhune]

I did a test of AddGeometryAttributes_management() for "POINT_X_Y_Z_M" and it can take 15 minutes to do an SDE feature class with about 80,000 points. I'm not even trying anything fancy with reprojecting, I'm just using the feature class's assigned coordinate system. The output is just Point_X and Point_Y fields. Is this typical performance?

[XanderBakker]

Sounds slow... Will do a test on a local dataset that I have and see how it compares to using the da update cursor.