""" Spatial Data Modeller for ESRI* ArcGIS 9.2 Copyright 2007 Gary L Raines, Reno, NV, USA: production and certification Don L Sawatzky, Spokane, WA, USA: Python software development """ import sys, string, os, traceback, math from floatingrasterarray import FloatRasterSearchcursor # Create the Geoprocessor object import arcgisscripting gp = arcgisscripting.create() # Check out any necessary licenses gp.CheckOutExtension("spatial") # Load required toolboxes... #gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx") gp.OverwriteOutput = 1 # Local variables... #Parameters PostProb = gp.GetParameterAsText(0) PPStd = gp.GetParameterAsText(1) TrainSites = gp.GetParameterAsText(2) UnitArea = gp.GetParameter(3) SaveFile = gp.GetParameterAsText(4) #Functions def ZtoF(Z): ##Function ZToF(Z As Double) As Double ##' SUBROUTINE ZTOF(Z, PZ) ##'c ##'C ... SUBROUTINE TO COMPUTE FREQUENCY FROM Z. EQ. 26.2.17 IN ##'c ##'C ABRAMOWITZ, M. AND STEGUN, I.A. (EDITORS) ##'c "HANDBOOK OF MATHEMATICAL FUNCTIONS" ##'C WITH FORMULAS, GRAPHS, AND MATHEMATICAL TABLES" ##'C PUB. BY NATIONAL BUREAU OF STANDARDS OF THE ##'C U.S. DEPARTMENT OF COMMERCE, 1964. ##'c ##'C CALLED BY SUBROUTINES COMP, REORD, & WDIST ##'c ##' DATA PI/3.141592654/, CONST6/0.2316419/, B1/0.319381530/, ##' + B2/-0.356563782/, B3/1.781477937/, B4/-1.821255978/, ##' + B5/1.330274429/ PI = 3.141592654 CONST6 = 0.2316419 B1 = 0.31938153 B2 = -0.356563782 B3 = 1.781477937 B4 = -1.821255978 B5 = 1.330274429 ##'c ##' x = Z Z = float(Z) X = Z ##' IF (Z.LT.0.0) X = -Z if X < 0.0 : X = -Z ##' t = 1# / (CONST6 * x + 1#) t = 1.0 / (CONST6 * X + 1.0) ##' PID = 2# * PI PID = 2.0 * PI ##' XX = -x * x / 2# XX = -X * X / 2.0 ##' XX = Exp(XX) / SQRT(PID) XX = math.exp(XX) / math.sqrt(PID) ##' PZ = (B1 * T) + (B2 * T * T) + (B3 * T**3) + (B4 * T**4) + ##' + (B5 * T**5) PZ = (B1 * t) + (B2 * t * t) + (B3 * (t ** 3)) + (B4 * (t ** 4)) + (B5 * (t ** 5)) ##' PZ = 1# - (PZ * XX) PZ = 1.0 - (PZ * XX) ##' IF (Z.LT.0.0) PZ = 1.0 - PZ if Z < 0 : PZ = 1.0 - PZ return PZ ##' Return ##' End try: basename = os.path.basename(PostProb) sdmuc = basename.split("_")[0] #gp.AddMessage(sdmuc) #CellSize CellSize = float(gp.CellSize) ExpNumTP = gp.GetCount(TrainSites) #Num of Selected sites #gp.AddMessage("%s = %s" % ("ExpNumTP",ExpNumTP)) ConvFac = (CellSize**2)/1000000.0 / UnitArea #gp.AddMessage("%s = %s" % ("ConvFac",ConvFac)) #Single values #PredT = SUM(Tras.Value) #gp.MakeRasterLayer_management(PostProb,"PostProbRL") ## desc = gp.describe(PostProb) ## if desc.datatype == 'RasterLayer': PostProb =desc.catalogpath TRasRows = FloatRasterSearchcursor(gp,PostProb) PredT = 0.0 for TRasRow in TRasRows: #gp.AddWarning("%s * %s = %s" % (TRasRow.Value,TRasRow.Count,TRasRow.Value*TRasRow.Count)) PredT += (TRasRow.Value * TRasRow.Count) #gp.AddMessage("%s = %s" % ("Sum N*P",PredT)) PredT *= ConvFac #gp.AddMessage("%s = %s" % ("Predicted T",PredT)) del TRasRow,TRasRows ## desc = gp.describe(PPStd) ## if desc.datatype == 'RasterLayer': PPStd =desc.catalogpath TRasRows = FloatRasterSearchcursor(gp,PPStd) TVar = 0.0 for TRasRow in TRasRows: #gp.AddWarning("%s * %s = %s" % (TRasRow.Value,TRasRow.Count,TRasRow.Value*TRasRow.Count)) TVar += (TRasRow.Value * TRasRow.Count * ConvFac )**2 #gp.AddMessage("%s = %s" % ("TVar",TVar)) TStd = math.sqrt(TVar) #gp.AddMessage("%s = %s" % ("TStd",TStd)) del TRasRow,TRasRows #TS = (PredT - ExpNumTP)/TStd TS = (PredT - ExpNumTP) / TStd #gp.AddMessage("%s = %s" % ("TS",TS)) #PostProb n = ExpNumTP T = PredT #STD = TStd P = ZtoF(TS) *100.0 if P>= 50.0: overallCI = 100.0 * (100.0 - P) / 50.0 else: overallCI = 100.0 * (100.0 - (50 + (50 - P))) / 50.0 Text = """ Overall CI: %(0).1f%%\r Conditional Independence Test: %(1)s\r Observed No. training pts, n = %(2)i\r Expected No. of training points, T = %(3).1f\r Difference, T-n = %(4).1f\r Standard Deviation of T = %(5).3f\r \r ------------------------------------------------\r Conditional Independence Ratio: %(6).2f \r Values below 1.00 may indicate conditional dependence\r among two or more of your data sets. \r \r ------------------------------------------------\r Agterberg & Cheng Conditional Independence Test\r \r This is a one-tailed test of the null hypothesis that T-n=0. The test\r statistic is (T-n)/standard deviation of T. Probability values greater\r than 95%% or 99%% indicate that the hypothesis of CI should be rejected,\r but any value greater than 50%% indicates that some conditional\r dependence occurs>\r \r Probability that this model is not conditionally independent with\r (T-n)/Tstd = %(7).6f is %(8).1f%%\r ------------------------------------------------\r \r Input Data:\r Post Probability: %(9)s\r Post Probability Std Deviation: %(10)s\r Training Sites: %(11)s \f """ % {'0': overallCI, '1': sdmuc, '2':n, '3':T, '4':T-n, '5':TStd, '6':n/T, '7':TS, '8':ZtoF(TS)*100.0, '9':PostProb, '10':PPStd, '11':TrainSites} #gp.AddMessage('Overall CI: %.1f%%'%overallCI) gp.AddMessage(Text) #PPprefix = os.path.basename(PostProb) if SaveFile: #savefile = os.path.join(gp.Workspace,'ACCITest_%s.txt'%SaveFile) file = open(SaveFile,"w") file.write(Text) gp.AddMessage("Text File saved: %s"%SaveFile) except: # get the traceback object tb = sys.exc_info()[2] # tbinfo contains the line number that the code failed on and the code from that line tbinfo = traceback.format_tb(tb)[0] # concatenate information together concerning the error into a message string pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n " + \ str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n" # generate a message string for any geoprocessing tool errors msgs = "GP ERRORS:\n" + gp.GetMessages(2) + "\n" gp.AddError(msgs) # return gp messages for use with a script tool gp.AddError(pymsg) # print messages for use in Python/PythonWin print pymsg print msgs raise