You can define an output projection (and transform) for a cursor and write out projected coordinates.# 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)
