Use arcpy to calculate areas in a different coordinate system

[AnthonyCheesman1]

Hi

I have an arcpy process generates a summarised dataset for me from a source dataset. Due to the limitations of our corporate spatial library, the source data (and my summarised data) are geographic (GDA94).

My summarised dataset has an AREA_HA field (area, hectares) and I would like to automatically populate that field as part of the summary process. I would like this calculation to be performed in the GDA94 VG94 projected coordinate system (shape@hectares doesn't provide the accuracy I need).

I'm a little at a loss as to how to do this programatically. I've had a look at spatial references (arcpy.SpatialReference) but it doesn't seem to do what I am looking for.

Can anyone advise how I would go about this?

Thanks

[DanPatterson_Retired · Accepted solution]

As long as you know what coordinate system you want the values in, you can use the existing tools in arcmap in the Add Geometry Attributes—Help | ArcGIS for Desktop  tool, the code example will show that you can use an optional spatial reference rather than the coordinates of the existing data

[AnthonyCheesman1]

Thanks Dan, I never knew that toolset existed.

Is there a way to capture the area value and write it directly to a pre-existing field, rather than have it create a new field?

The other option would be create the new field, copy the new field to the existing field, and then delete the new field I suppose.

[DanPatterson_Retired]

Yes...via code...my suggestions would be to leave fields that have calculated values intact, and produce the results of a new operation in a new field thereby preserving the legacy information.  Once could also string this tool with delete field in a small scriptlet should one really really be adverse to having both sets of data there.  It is a matter of preference and desired work expenditure.

[RebeccaStrauch__GISP]

As Dan mentioned, it's sometimes nice to keep track of the x/y from both coordinaste systems.  In the project I am working on right now, I need the x/y in my projected system, and in DD, DDM, and DMS, and I'm keeping two feature classes for the field.  Below I've included a simplified snippet of my script (which selects polygons, creates random points, grabs elevation, and added all the coordinates...more than is needed here).  This won't run as is, since I have some hardcoded variab les below that you will have to work with, but you should be able to see how you can go about adding the additional fields.

Hope this helps.  (....EDIT: added code

# Import modules
import time
import arcpy
from arcpy import env

# Set the Geoprocessing environment...Hardcoded to NAD83....change if needed
geoSR = arcpy.SpatialReference(4269)

# Script arguments...
ProjDir = r"C:\_beartest"
studyPoly = r"C:\_beartest\Prep.gdb\Unit20A_edit"
stydGDB,stdyFile = os.path.split(studyPoly)
arcpy.env.workspace = stydGDB
arcpy.env.overwriteOutput = True 

# Local variables...
#   grabs the Spatial Reference of the study area, so save to fiel, and for reproject
srDesc = arcpy.Describe(studyPoly).spatialReference
srName = "'{0}'".format(srDesc.name)      #format ends up being "'the_text_here'"  notice "' '"
srWkid = "'{}'".format(srDesc.PCSCode)
print("  ->Initial Projection of {0}: {1}; WKID: {2}".format(stdyFile, srName, srWkid))

#  prep.gdb is typically where this is written, with the project folder...
prepFGDB =   ("{0}\prep.gdb".format(ProjDir))

#  there are many temp FCs that are created in this process
randomPtSet     = arcpy.os.path.join(prepFGDB, "Pts{0}b".format(TransSeriesNo))  # will add elevation value to this
finaltmp1X      = arcpy.os.path.join(prepFGDB, "Pts{0}1X".format(TransSeriesNo))
finaltmp2X      = arcpy.os.path.join(prepFGDB, "Pts{0}2X".format(TransSeriesNo))    
#
finalOutdd      = arcpy.os.path.join(prepFGDB, "Pts{0}_dd".format(TransSeriesNo))
finalOutProj    = arcpy.os.path.join(prepFGDB, "Pts{0}_{1}".format(TransSeriesNo, srWkid.strip("'")))
#
Point_Count = "0"     # so point number for the series starts at Series number * 100

arcpy.AddXY_management(randomPtSet)  

# create list of fields to add..    
lstAddFields = [["projtd_x", "DOUBLE", "13", "5", "", "x"], 
                  ["projtd_y", "DOUBLE", "13", "5", "", "y"], 
                  ["xyProj", "TEXT", "", "", "50", "xyProjection"],
                  ["xyProjID", "SHORT", "6", "", "", "xyWKID"],  
                  ["ptdd_lat", "DOUBLE", "9", "7", "", "lat_dd"], 
                  ["ptdd_long", "DOUBLE", "11", "7", "", "long_dd"]] 

# Add fields from list above
print("->Adding Fields....")
for addFld in lstAddFields:
    print("   ...{}".format(addFld[0]))
    arcpy.AddField_management(randomPtSet, addFld[0], addFld[1], addFld[2], addFld[3], addFld[4], addFld[5], "#", "NON_REQUIRED")


# grabbing X/Y coordinates and the source study area projection name and wkid
arcpy.CalculateField_management(randomPtSet, "projtd_x", "!POINT_X!", "PYTHON_9.3", "")
arcpy.CalculateField_management(randomPtSet, "projtd_Y", "!POINT_Y!", "PYTHON_9.3", "")
arcpy.CalculateField_management(randomPtSet, "xyProj", srName, "PYTHON_9.3", "")
arcpy.CalculateField_management(randomPtSet, "xyProjID", srWkid, "PYTHON_9.3", "")

# adding POINT_X and POINT_Y geographic (long/lat) coords WITHOUT projecting...
#    and calculate these DD fields so Long has negative sign instead of W at end
arcpy.AddGeometryAttributes_management(randomPtSet, "POINT_X_Y_Z_M", "", "", geoSR) # "4269")
arcpy.CalculateField_management(randomPtSet, "ptdd_long", "!POINT_X!", "PYTHON_9.3", "")
arcpy.CalculateField_management(randomPtSet, "ptdd_lat", "!POINT_Y!", "PYTHON_9.3", "")

# NOTE: ConvertCoordinateNotation actually PROJECTs the source and outputs a new layer
#  http://desktop.arcgis.com/en/desktop/latest/tools/data-management-toolbox/convert-coordinate-notatio...
# The tool projects the output to the spatial reference specified. If the input and output coordinate 
#  systems are in a different datum, then a default transformation is used based on the coordinate 
#  systems of the input and the output and the extent of the data.

# Converting the random points to the geoSR hardcoded around line 89
# NOTEL commented line would add in same format as the calc field above, but same default fieldnames as DD_2, so
#   would be DDLat_1 and DDLong_1...I prefer to set to ptdd long and ptdd_lat rather than have the _1
#      arcpy.ConvertCoordinateNotation_management(finaltmp1X, finaltmp2X, "POINT_X", "POINT_Y", "DD_2", "DD_NUMERIC", "", geoSR )  
# So calculate:
#    to decimal degrees with N/S/E/W designation at end (vs. ptdd_ versions above that have W/S as negative)
arcpy.ConvertCoordinateNotation_management(randomPtSet, finaltmp1X, "POINT_X", "POINT_Y", "DD_2", "DD_2", "", geoSR )  
#    to degree-decimal minutes
arcpy.ConvertCoordinateNotation_management(finaltmp1X, finaltmp2X, "POINT_X", "POINT_Y", "DD_2", "DDM_2", "", geoSR )  
#    to degree-minutes-seconds
arcpy.ConvertCoordinateNotation_management(finaltmp2X, finalOutdd, "POINT_X", "POINT_Y", "DD_2", "DMS_2", "", geoSR) 

# projecting back to initial, so have projects and DD versions
print("finalOutProj: {0}".format(finalOutProj))
arcpy.Project_management(finalOutdd, finalOutProj, srDesc) 

# Process: Delete temp files...
tempFiles = [randomPtSet, finaltmp1X, finaltmp2X]
for tmp in tempFiles:
    if arcpy.Exists(tmp):
        arcpy.Delete_management(tmp)

print(' \n !!! Success !!!  ')

[RebeccaStrauch__GISP]

As an addendum to my code above, I've moved on to my next step in which case I grab the coordinates and other fields.  In testing, I'm noticing the "SHAPE@XY" of both the projected and geographic point featureclass do not match up to the stored values.  I'm looking into it, and my guess is it is a projection issue.  Not exactly sure why (using standard wkid) but will put in a tech support ticket.  just fyi.