Geoprocessing
From DreamsteepWiki
I wasted a year of my life learning this ESRI crap, do yourself a favor and go learn Opensource tech like postGIS instead!
READ THIS ----> Setup_database_server_on_ubuntu IT WILL CHANGE YOUR LIFE!
IF you insist on going down that path check out my:
geometry processing arcgeo table processing arcattr spatial processing spatial_treatment
GEOPROCESSING NOTES
arcgisscripting module
arcgisscripting is a python bytecode module , add it to your PYTHONPATH it is in the esri/arcgis/bin folder
instantiate the geoprocessor ( make sure win32.com in installed)
import win32com.client
gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")
BATCH
#merge shapefiles
import arcgisscripting
import os
files = ['C:/GIS_DATA/ee/test.shp','C:/GIS_DATA/grid_park_co.shp']
gp = arcgisscripting.create(9.3)
outpath = 'C:/GIS_DATA'
###
gp.Workspace = outpath
gp.toolbox = "management"
# Set a variable to store the updated feature class
out_feat_class = "mass_towns.shp"
if os.path.lexists((outpath+'/'+out_feat_class) )==0:
gp.CreateFeatureclass( outpath , out_feat_class ,'POLYGON' ) #,spatref)
for file in files:
if os.path.lexists(file):
gp.Append( file , (outpath+'/'+out_feat_class),"NO_TEST" )
if os.path.lexists(file)==0:
print file + ' DOES NOT EXIST '
GEODATABASE
A little rant:
Geodatabase is a SCAM!
ESRI is trying to snare people into its own custom formats,GEODATABASES ARE CRAP! ESRI is notorious for breaking a single program into "modules" and then nickel and diming its customers to death (mostly education and government who aren't spending their own money so nobody complains) . Mr. Jack Dangermond , take your gp.checkOutExtension and shove it up your greedy arse!
create and batch populate
# Import native arcgisscripting module
#
import arcgisscripting
import sys
import os
# Create the geoprocessor object
gp = arcgisscripting.create(9.3)
# Allow for the overwriting of file geodatabases, if they previously exist
gp.OverWriteOutput = 1
# Set workspace to folder containing personal geodatabases
#gp.Workspace = sys.argv[1]
# Identify personal geodatabases
#pgdbs = gp.ListWorkspaces()
pgdb = sys.argv[1]
# Set workspace to current personal geodatabase
gp.workspace = pgdb
# Create file geodatabase based on personal geodatabase
fgdb = pgdb[:-4] + ".gdb"
if os.path.lexists(os.path.dirname(fgdb))==0:
gp.CreateFileGDB(os.path.dirname(fgdb), os.path.basename(fgdb))
# Identify feature classes and copy to file gdb
fcs = gp.ListFeatureClasses()
# Identify tables and copy to file gdb
tables = gp.ListTables()
# Identify datasets and copy to file gdb
# Copy will include contents of datasets
datasets = gp.ListDatasets()
print '###########'
print os.path.basename(fgdb)
print fcs
print tables
print datasets
raw_input()
for fc in fcs:
print "Copying feature class " + fc + " to " + fgdb
gp.Copy(fc, os.path.join(fgdb, fc))
#for table in tables:
# print "Copying table " + table + " to " + fgdb
# gp.Copy(table, os.path.join(fgdb, table))
for dataset in datasets:
print "Copying dataset " + dataset + " to " + fgdb
gp.Copy(dataset, os.path.join(fgdb, dataset))
query type
# Create search cursor
rows = GP.SearchCursor("C:/data/franky.shp")
row = rows.Next()
# Calculate the total length of all roads
while row:
# Create the geometry object
feat = row.shape
shapetype= str(feat.Type)
#
if shapetype == 'polygon':
print 'POLYGON DETECTED \n'
#
if shapetype == 'polyline':
print 'LINE DETECTED \n'
#
if shapetype == 'point':
print 'POINT DETECTED \n'
#length = length + feat.Length
row = rows.Next()
edit feature class
# Open an insert cursor for the feature class
cur = GP.InsertCursor(fcname)
# Create the array and point objects needed to create a feature
lineArray = GP.CreateObject("Array")
pnt = GP.CreateObject("Point")
# Add two points to the array
pnt.x = 358331
pnt.y = 5273193
lineArray.Add(pnt)
pnt.x = 358337
pnt.y = 5272830
lineArray.Add(pnt)
# Create a new row, or feature, in the feature class
feat = cur.NewRow()
# Set the geometry of the new feature to the array of points
feat.shape = lineArray
# Insert the feature
cur.InsertRow(feat)
error handling
try:
#Input Data
fromFC = r"C:/data/keith.mdb/may"
except:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
pymsg = "PYTHON ERRORS:\nTraceback Info:\n%s\nError Info:\n %s: %s\n" % \
((tbinfo,) + sys.exc_info()[:2])
AddPrintMessage(pymsg, 2)
msgs = "GP ERRORS:\n" + gp.GetMessages(2) + "\n"
AddPrintMessage(msgs, 2)
del tb
standard esri messaging function
def AddPrintMessage(msg, severity):
print msg
if severity == 0: gp.AddMessage(msg)
elif severity == 1: gp.AddWarning(msg)
elif severity == 2: gp.AddError(msg)
preform a desribe
desc = gp.describe(infc)
featuretype = desc.FeatureType
datatype = desc.DataType
print '# # DATATYPE IS ' + datatype + '\n'
print '# # FEATURTYPE IS ' + featuretype + '\n'
using a search cursor
from win32com.client import Dispatch
GP = Dispatch("esriGeoprocessing.GpDispatch.1")
# Create search cursor
rows = GP.SearchCursor(sys.argv[1])
row = rows.Next()
while row:
row = rows.Next()
load a formated text file and create a shapefile
"""
#EXAMPLE DATA
1;22.4;44.2 ;
1;55.4;56.3 ;
2;65.4;66.3 ;
2;75.4;76.3 ;
3;175.4;76.3 ;
3;275.4;-576.3 ;
###
1 ;422.4 ;44.2 ;
1 ;55.4 ;756.3 ;
2 ;565.4 ;66.3 ;
2 ;75.4 ;5476.3 ;
3 ;6175.4 ;776.3 ;
3 ;275.4 ;-576.3 ;
"""
import win32com.client ,sys,os
import sys_file_io #fileinput
gp = win32com.client.Dispatch('esriGeoprocessing.GpDispatch.1')
gp.workspace = 'C:/data'
gp.overwriteoutput = 1
def readFeatureFile ():
gp.CreateFeatureclass('C:/data','foo2.shp','POLYLINE')
#gp.InsertCursor()
cur = gp.InsertCursor("C:/data/foo2.shp")
lineArray = gp.CreateObject("Array")
pnt = gp.CreateObject("Point")
ID = -1
#for line in fileinput.input("C:/crap.txt"):
# values = line.split(";")
fhandler = sys_file_io.file_io_three()
fhandler.readfilelines('C:/crap.txt')
for line in fhandler.filecontents_list:
values = line.split(';')
#print values[0]
#print values[1]
#print values[2]
temp1 = values[0]
temp2 = values[1]
temp3 = values[2]
#print values
if len(temp1):
pnt.id = str(temp1)
pnt.x = str(temp2)
pnt.y = str(temp3)
#print pnt.id
#print pnt.x
#print pnt.y
lineArray.add(pnt)
if ID == -1:
ID = pnt.id
if ID!=pnt.id:
#new row or feature in feature class
feat = cur.NewRow()
#set geometry to array of points
feat.shape = lineArray
#insert feature
cur.InsertRow(feat)
lineArray.RemoveAll()
lineArray.add(pnt)
ID=pnt.id
readFeatureFile()
scan feature class row and field data
gp.Workspace = "C:/data"
fcs = gp.ListFeatureClasses("*")
#fcs = gp.ListRasters("*") #RasterDataset
#fcs = gp.ListDataSets('*')
fcs.reset()
fc = fcs.next()
#go through each feature class in the workspace
while fc :
desc = gp.Describe(fc)
#GET INFO ABOUT FEATURE CLASS
typeoffeature= str(desc.DataType )
#print fc +'\n'
fields = gp.ListFields(fc)
fields.reset()
field = fields.next()
print '\n'+ '#OBJECT IS '+ fc + '######################'
while field:
print '#FIELD NAME ' +field.Name
#print 'LENGTH ' +str(field.Length)
#print '\n'
field = fields.next()
#GET VALUES OF EACH FIELD IN EACH FEATURE CLASS
#find each field and print the value stored in it
if typeoffeature == 'RasterDataset':
print '#DEBUG# RASTER FILE DETECTED \n'
if typeoffeature == 'Tin':
print '#DEBUG# RASTER FILE DETECTED \n'
if typeoffeature == 'ShapeFile':
print '#DEBUG# SHAPE FILE DETECTED \n'
#ONLY ATTEMPT TO READ TABLES ON A SHAPE FEATURE (OTHERS WILL WORK ALSO)
if typeoffeature == 'ShapeFile' :
rows = gp.SearchCursor(fc)
row = rows.Next()
while row:
print '# #FID FIELD (EACH ROW) ##########\n'
print row.FID #DIRECTLY ACCESS A FIELD
#print row.GetValue('FID') #BY NAME
#print row.GetValue(field.Name)
row = rows.Next()
del row
del rows
##
fc=fcs.next()
read a feature as txt data
import win32com.client ,sys,os,string
import sys_file_io #fileinput
gp = win32com.client.Dispatch('esriGeoprocessing.GpDispatch.1')
gp.workspace = 'C:/data'
gp.overwriteoutput = 1
output = []
outfilename = 'C:/crap.txt'
infilename = 'C:/data/sixseven/67th.shp'
def readGeoTXT ():
rows = gp.SearchCursor( infilename )
row=rows.Next()
a=0
while row:
feat=row.shape
print '#### TYPE '+str(feat.Type) +'\n'
#print 'PART IS ' + '\n'
print 'PART COUNT IS '+str(feat.PartCount)
while a <feat.PartCount:
#get each part of geo
roadArray=feat.GetPart(a)
roadArray.Reset() #no () in example #?
#Get the first point
pnt = roadArray.Next()
while pnt:
#print str(pnt.id)+';'+str(pnt.x)+';'+str(pnt.y)
output.append( str(pnt.id)+';'+str(pnt.x)+';'+str(pnt.y)+' \n')
pnt= roadArray.Next()
a =a+1
row=rows.Next()
print 'CREATING FILE '+outfilename
fhandler = sys_file_io.file_io_three()
#fhandler.readfilelines('C:/crap.txt')
fhandler.writefile_list(outfilename,output)
## ## ##
readGeoTXT()
list feature fields
try:
#Get Inputs (feature class and feature id)
infc = gp.GetParameterAsText(0)
#poly_id = gp.GetParameterAsText(1)
print infc
#Describe input fc. Describe object contains useful information used later.
desc = gp.describe(infc)
featuretype = desc.FeatureType
datatype = desc.DataType
print '# # DATATYPE IS ' + datatype + '\n'
print '# # FEATURTYPE IS ' + featuretype + '\n'
except:
print 'command got hosed\n'
print 'hit return to exit\n'
raw_input()
