I have a script Tool developed by one of my previous colleague. It was running in a remote Computer so that everybody can use These Tools whenever needed. I was working well but suddenly it has stopped working and its a Panic. We have to fix the error of script tool immediately as we have only one 3d Analyst license in arcgis 9.3.
our python code is as follows:
try:
from osgeo import gdal
except:
sys.exit("No GDAL library!")
import string
import sys, os
import csv
import arcgisscripting
gp = arcgisscripting.create(9.3)
#overwriting True
gp.OverWriteOutput = 1
#--------------------------------------------------------------------------
# INPUTS
## test data
##i= "F:/Projekte/A1101003/Grundlagen/DSM/griend_dtm0_GK3.tif"
##shp_mask = "F:/Projekte/A1101003/ArcGis/Themen/Polderung_Var02/Block45_Polderung_Var02_DHDN.shp"
##rasterpath = "E:/test/pv3"
##ordinate = "Ordinate"
##fieldName = "Polder"
##cellsize= "0,25"
i= gp.GetParameterAsText(0) #raster
shp_mask = gp.GetParameterAsText(1) #shp with multiple features overlaping raster
fieldName = gp.GetParameterAsText(2) #field with name used for cliped rasters
ordinate = gp.GetParameterAsText(3) #field with height value as 'Ordinate' to create rasters
rasterpath = gp.GetParameterAsText(4) #output folder
filename=rasterpath+"/_CutFill.csv"
cellsize =gp.Describe(i).MeanCellHeight
if not gp.Describe(i).SpatialReference.Name==gp.Describe(shp_mask).SpatialReference.Name:
gp.AddError("Shapefile coordinate system is different from raster coordinate system!!!\
Please reproject the shapefile to raster coordinate system '%s'." %gp.Describe(i).SpatialReference.Name)
print gp.GetMessages()
raise Exception()
#--------------------------------------------------------------------------
# DEFINITIONS of FUNCTIONS
def createDic (shp, field_name, field_ordinate):
shpdict={}
rows=gp.SearchCursor(shp)
row=rows.Next()
while row:
#shapefile field name - the ordinate value
shpdict[str(row.GetValue(field_name))] = "%.2f" %row.GetValue(field_ordinate)
row=rows.Next()
del row, rows
return shpdict
def cutfill4clip (raster, shape, field_name, field_ordinate):
"""
Cut/Fill raster created from clipedDEM (every feat) and raster of Ordinate(every feat)
"""
gp.CheckOutExtension("3d")
gp.workspace= rasterpath
rows=gp.SearchCursor(shape)
row=rows.Next()
count = 1
global increment
while row:
#--- setting progressor
gp.SetProgressorPosition(int(round(count*increment)))
#--- setting names
# checking if name has a space, changing to '-'
row_name = str(row.GetValue(field_name))
if " " in row_name:
row_name = row_name.replace(" ", "-")
#output of cliped raster
output = "c_"+ row_name[:10]
#---------------------
#getting rounded value of Height
if isinstance(row.GetValue(field_ordinate), float):
value=row.GetValue(field_ordinate) #float value ex:5.520000000001
rvalue= "%.2f" %value #getting string as output
rvalue= rvalue.replace(".", ",")
else:
value=row.GetValue(field_ordinate).replace(".", ",")
rvalue = value
#-----------------------
#output of ordinate raster
output_ordinate = "ord_"+str(row_name)[:9]
#output of cut/fill raster
output_cutfill = "cf_"+str(row_name)[:9]
#RasterClip_management
clipgeo=row.Shape
clipextent=str(row.Shape.Extent)
try:
gp.clip_management(raster, clipextent, output, clipgeo ,"0", "ClippingGeometry")
except:
msg = gp.GetMessages()
gp.AddMessage(msg)
sys.exit()
gp.AddMessage("Raster clipped")
#CalculateStatistics_management(in_raster_dataset, x_skip_factor, y_skip_factor, ignore_values)
gp.CalculateStatistics_management(output, "1", "1", "0")
#Creating an raster of height from ordinate
extent= clipextent[:-16] #no NaN NaN NaN NaN values
dist = "uniform "+rvalue+" "+rvalue #special format for raster creation
gp.CreateRandomRaster_management(rasterpath, output_ordinate, dist, extent, cellsize)
#Creating cutfill raster
gp.CutFill_3d(output, output_ordinate, output_cutfill, '')
gp.AddMessage("... cut/fill raster %s" %output)
row=rows.Next()
count+= 1
del row, rows
#--------------------------------------------------------------------------
# MAIN
#
#--- setting Progressor
result = gp.GetCount_Management(shp_mask)
num_of_rows = int(result.GetOutput(0))
increment = 100.0/num_of_rows
gp.SetProgressor("step", "Clipping rasters with Shp ...", 0, 100, 1)
gp.AddMessage("Clipping raster with Shapefile and creating Cut/Fill rasters.")
if os.path.exists(shp_mask):
cutfill4clip (i, shp_mask, fieldName, ordinate) #goes to rasterpath (defined by user)
else:
print "Shapefile inputs are invalid"
#...creating a dictionary for csv
rdict = createDic (shp_mask, fieldName, ordinate)
#...writing to CSV
fCal=open(filename, 'wb') #binary opened!!!
csv.register_dialect('dotcoma', delimiter=';') #needed to open with Excel
writer=csv.writer(fCal, dialect='dotcoma')
#writing first row
writer.writerow(('CF_raster', ''+fieldName+'' , 'Ordinate', 'Volume_Pos', 'Volume_Neg', 'Area_Pos'))
gp.AddMessage("Going through rasters to create a csv table.")
gp.workspace = rasterpath
gp.CheckOutExtension ("3D")
#...looping through rasters of ('cf_...') to fill CSV
for raster in gp.listdatasets("cf_*"):
rasterName = raster.split("_", 1)[1].encode()
if rasterName not in rdict:
rasterName = rasterName.capitalize()
if rasterName not in rdict:
rasterName = rasterName.replace("-", " ")
if rasterName not in rdict:
gp.AddMessage('Shit!')
gp.SetProgressor("default", "getting data for Polders...")
gp.AddMessage('... getting data for Polder: %s' %rasterName)
rows=gp.SearchCursor(raster)
row=rows.Next()
vol_pos, vol_neg = 0, 0
while row:
if row.VOLUME < 0: # was area- 200 m2
vol_neg+=row.VOLUME
if row.VOLUME > 0:
vol_pos+=row.VOLUME
row=rows.Next()
#preparing for CSV (. changing to ,)
vol_neg, vol_pos = str(vol_neg).replace(".", ","), str(vol_pos).replace(".", ",")
ord_val = str(rdict[rasterName]).replace(".", ",")
#getting the Area 3D Value above the Ordinate
gp.toolbox = "3D"
gp.surfacevolume_3d ("c_"+raster[3:], "#", "ABOVE", ord_val, "#")
msg=gp.GetMessages() # could write somewhere...
#msg conatains: 2D, 3D, Volume
#searching for 2d volume
start = (string.find(msg, "2D Area="))+8
end = string.find(msg, "3D", start)
area_pos = msg[start:end].strip()
#writing just for every raster all PositiveVolume and all NegativeVolume
writer.writerow((raster, rasterName, ord_val, vol_pos, vol_neg, area_pos))
del row, rows
fCal.close()
#--------------------------------------------------------------------------
os.startfile (filename)
**The process Shows a Progress of 6% and then stop working and
the error Message printed during the processing is like this:**
Executing: CutFillDEMOrdinate F:\Projekte\A1001025\ArcGis\Themen-Planung\B-Hoehe\hoehe_hn+50 F:\Projekte\A1001025\ArcGis\Themen-Planung\Polderplanung_Var6_Merc.shp FID Ordinate F:\Projekte\A1001025\ArcGis\Themen-Planung\temp
Start Time: Wed Apr 02 09:35:53 2014
Running script CutFillDEMOrdinate...
Clipping raster with Shapefile and creating Cut/Fill rasters.
Executing: Clip F:\Projekte\A1001025\ArcGis\Themen-Planung\B-Hoehe\hoehe_hn+50 "3413253,28963566 5878810,1020835 3413413,96707874 5878983,76715465" F:\Projekte\A1001025\ArcGis\Themen-Planung\temp\c_0 in_memory\fC0D8CD67_0B59_4BBE_ACDD_62DD737B7EDF 0 ClippingGeometry
Start Time: Wed Apr 02 09:35:56 2014
ERROR 999998: Unexpected Error.
Failed to execute (Clip).
End Time: Wed Apr 02 09:35:59 2014 (Elapsed Time: 3,00 seconds)
<type 'exceptions.SystemExit'>:
Failed to execute (CutFillDEMOrdinate).
End Time: Wed Apr 02 09:35:59 2014 (Elapsed Time: 6,00 seconds)
**Last few days we are trying to solve it. I did not use also arcgis 9.3 Long time. If someone could help or recomend a possible solution will be a great help for us**
Thanks a lot
Shiuli