Use arcpy to calculate areas in a different coordinate system

2703
5
Jump to solution
11-03-2015 04:35 PM
AnthonyCheesman1
Occasional Contributor II

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

0 Kudos
1 Solution

Accepted Solutions
DanPatterson_Retired
MVP Esteemed Contributor

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

View solution in original post

5 Replies
DanPatterson_Retired
MVP Esteemed Contributor

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

View solution in original post

AnthonyCheesman1
Occasional Contributor II

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.

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

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
MVP Esteemed Contributor

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
MVP Esteemed Contributor

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.

0 Kudos