Reproject a point with arcpy?

DennisGeasan

I would like to use arcPy to reproject XY point fields in a table to a different coordinate system.  Essentially do the same as the Calculate Geometry option but from a Python script.  Is it possible?  I can't find any functions in arcPy that provides reprojection of the point values derived from the 'shape' field or from XY values stored as attributes.

DG

KimOllivier

There is a projection parameter option on the cursor objects that will reproject to new coordinates as you read the records. Getting a transform to work as well is a bit tricky, but basically set that as well as a string.

# extent_xy.py
# add extents of parcels in NZMG to a NZTM dataset for backward compatibility
# update to use on the fly projection directly from NZTM layers
# see Werner Flacke's book and script example
# updated for 9.2
# Kim Ollivier 13 April 2007
# input NZTM coverage library, parcel and plabel coverages
# output csv files to load and join back to coverages using AML
# issues solved:
#               add spatial reference to SearchCursor
#               define a Geographic Transformation that works with a coverage
#               because coverages cannot define NZGD2000 use WGS84 definition
#               or the transform is ignored

import os,sys,glob,arcgisscripting

gp = arcgisscripting.create()

def parcel(tile) :
    print "Parcel",tile
    ws = "e:/lib/nztm/tile/t"+str(tile)+"/data"
    outfolder = "e:/lib/nztm/tile/t"+str(tile)+"/data"
    ds = ws + "/parcel"
    f1 = open(outfolder+"/e_parcel.txt","w")
    f1.write("par_id,eminx_nzmg,eminy_nzmg,emaxx_nzmg,emaxy_nzmg\n")
    gp.Workspace = ds
    print ds
    # this is a good on-the-fly switch for the searchcursor
    # srOut must be a com object, not a file ref or a factory code
    rows = gp.SearchCursor(ds+"/polygon",'"PAR_ID" > 0',srOut) 
    n = 0
    rows.Reset()
    row = rows.Next()
    while row  :
        print >> f1,"%d,%s" % (row.par_id,row.shape.Extent.replace(" ",","))
        row = rows.Next()
        n +=1
    del rows
    print n,"Polygons"
    f1.close()
    return

def plabel(tile) :
    print "Plabel",tile
    ws = "e:/lib/nztm/tile/t"+str(tile)+"/data"
    ds = ws + "/plabel"
    outfolder = "e:/lib/nztm/tile/t"+str(tile)+"/data"
    f2 = open(outfolder+"/e_plabel.txt","w")
    f2.write("par_id,eminx_nzmg,eminy_nzmg\n")
    gp.Workspace = ds
    rows = gp.SearchCursor(ds+"/point",'"PAR_ID" > 0',srOut)
    n = 0
    rows.Reset()
    row = rows.Next()
    while row  :
        n += 1
        print >>f2,"%d,%s,%s" % (row.par_id,row.shape.Extent.split()[0],row.shape.Extent.split()[1])
        row = rows.Next()
    del rows
    print n,"plabels"
    f2.close()
    return

# ------ main ----
#
srOut = gp.CreateObject("SpatialReference")
srOut.CreateFromFile("c:/arcgis/nzmg.prj")
# WGS_1984 works for coverages defined with an equivalent custom Transverse projection definition
gp.GeographicTransformations = 'NZGD_1949_To_WGS_1984_3_NTv2' # ;New_Zealand_1949_To_NZGD_2000_3_NTv2'
# but is this being used?
inFC = "e:/lib/nztm/tile/t1009/data/plabel/point"
# not at 9.3?? outFC = "in_memory/dummy1.shp"
outFC = "c:/tmp/dummy1.shp"
outPrj = "c:/arcgis/nzmg.prj" # cannot be a COM object or factory id 43040
geoTrans = "NZGD_1949_To_WGS_1984_3_NTv2" # ;New_Zealand_1949_To_NZGD_2000_3_NTv2"
if gp.Exists("c:/tmp/dummy1.shp") :
    gp.Delete("c:/tmp/dummy1.shp")
#
gp.OverwriteOutput = 1
gp.MakeFeatureLayer_management(inFC,"tlayer","PLABEL# < 3")
gp.Project_management("tlayer",outFC,outPrj,geoTrans)

try :
    if sys.argv[1].upper() == 'ALL' :
        lstTile = range(1001,1013)
    else :
        lstTile = [ tile for t in sys.argv[1].split(",")]
except :
    lstTile = range(1001,1013)
print "begin extent_xy"    
for t in lstTile :
    print t
    parcel(t)
    plabel(t)