# Name: make_fast_mosaic_v93.py # # Description # ----------- # This script will mosaic a bunch of grids, but using a secret method that is # mucho faster than the standard v9.3.1 SP1 mosaic tool. # # Written By: Chris Snyder, WA DNR, 03/01/2010, chris.snyder(at)wadnr.gov # # Written For: Python 2.5.1 and ArcGIS v9.3.1 SP1 # # UPDATES: # # Notes on input parameters (for the toolbox): # VARIABLE PAREMETER_INDEX PARAMETER_DATA_TYPE # ------------------------------------------------------------------- # gridList 0 Grids # outputGrid 1 Grid # tempWorkspace 2 Workspace # mosaicMethod 3 String (MAX, MIN, or MEAN) # tempNoDataValue 4 Long Integer try: #Process: Import some modules import os, string, sys, time, traceback, arcgisscripting #Process: Create the gp object gp = arcgisscripting.create(9.3) #Process: Defines some functions used for getting messages from the gp and python def showGpMessage(): gp.AddMessage(gp.GetMessages()) print >> open(logFile, 'a'), gp.GetMessages() print gp.GetMessages() def showGpWarning(): gp.AddWarning(gp.GetMessages()) print >> open(logFile, 'a'), gp.GetMessages() print gp.GetMessages() def showGpError(): gp.AddError(gp.GetMessages()) print >> open(logFile, 'a'), gp.GetMessages() print gp.GetMessages() def showPyLog(): #just print to the log file! print >> open(logFile, 'a'), str(time.ctime()) + " - " + message def showPyMessage(): gp.AddMessage(str(time.ctime()) + " - " + message) print >> open(logFile, 'a'), str(time.ctime()) + " - " + message print str(time.ctime()) + " - " + message def showPyWarning(): gp.AddWarning(str(time.ctime()) + " - " + message) print >> open(logFile, 'a'), str(time.ctime()) + " - " + message print str(time.ctime()) + " - " + message def showPyError(): gp.AddError(str(time.ctime()) + " - " + message) print >> open(logFile, 'a'), str(time.ctime()) + " - " + message print str(time.ctime()) + " - " + message #Specifies the root directory variable, defines the logFile variable, and does some minor error checking... dateTimeString = str(time.strftime('%Y%m%d%H%M%S')) scriptName = os.path.split(sys.argv[0])[-1].split(".")[0] userName = string.lower(os.environ.get("USERNAME")).replace(" ","_").replace(".","_") tempPathDir = os.environ["TEMP"] logFileDirectory = r"\\snarf\am\div_lm\ds\gis\tools\log_files" if os.path.exists(logFileDirectory) == True: logFile = os.path.join(logFileDirectory, scriptName + "_" + userName + "_" + dateTimeString + ".txt") try: print >> open(logFile, 'a'), "Write test successfull!" except: logFile = os.path.join(tempPathDir, scriptName + "_" + userName + "_" + dateTimeString + ".txt") else: logFile = os.path.join(tempPathDir, scriptName + "_" + userName + "_" + dateTimeString + ".txt") if os.path.exists(logFile)== True: os.remove(logFile) message = "Created log file " + logFile; showPyMessage() message = "Running " + sys.argv[0]; showPyMessage() #Process: Attempts to check out a Spatial Analyst license try: gp.CheckOutExtension("Spatial") except: message = "ERROR: Spatial analyst license is unavailable... Exiting script!"; showPyError(); sys.exit() #Process: Check out the highest license available try: if gp.CheckProduct("ArcView") == "Available": gp.SetProduct("ArcView") elif gp.CheckProduct("ArcEditor") == "Available": gp.SetProduct("ArcEditor") elif gp.CheckProduct("ArcInfo") == "Available": gp.SetProduct("ArcInfo") except: message = "ERROR: Could not select an ArcGIS license level! Exiting script..."; showPyError(); sys.exit() message = "Selected an " + gp.ProductInfo() + " license"; showPyMessage() #Process: Sets some gp environment variables gp.overwriteoutput = True gp.pyramid = "NONE" gp.rasterStatistics = "NONE" #Process: Collect the input parameters ## gridList = r"D:\csny490\usgs_dem10m\dem11746;D:\csny490\usgs_dem10m\dem11747" ## outputGrid = r"D:\csny490\temp\temp\howaboutthis" ## processingWorkspace = r"D:\csny490\temp\temp" ## mosaicMethod = "MAX" ## tempNoDataValue = "-32000" gridList = gp.GetParameterAsText(0) outputGrid = gp.GetParameterAsText(1) processingWorkspace = gp.GetParameterAsText(2) mosaicMethod = gp.GetParameterAsText(3) tempNoDataValue = gp.GetParameterAsText(4) #Process: Print out the input parameters message = "INPUT PARAMETERS"; showPyMessage() message = "----------------"; showPyMessage() message = "Input Grids = " + gridList; showPyMessage() message = "Output Grid = " + outputGrid; showPyMessage() message = "Processing Workspace = " + processingWorkspace; showPyMessage() message = "Mosaic Method = " + mosaicMethod; showPyMessage() message = "Temp NoData Value = " + tempNoDataValue + "\n"; showPyMessage() #Process: Figure out the maximum extent of all the input grids message = "Calculating maximum extent of input rasters..."; showPyMessage() i = 0 gridsToProcessList = gridList.split(";") for grid in gridsToProcessList: i = i + 1 dsc = gp.describe(grid) if i == 1: xMin = dsc.extent.xmin yMin = dsc.extent.ymin xMax = dsc.extent.xmax yMax = dsc.extent.ymax else: if dsc.extent.xmin < xMin: xMin = dsc.extent.xmin if dsc.extent.xmax > xMax: xMax = dsc.extent.xmax if dsc.extent.ymin < yMin: yMin = dsc.extent.ymin if dsc.extent.ymax > yMax: yMax = dsc.extent.ymax gp.extent = str(xMin) + " " + str(yMin) + " " + str(xMax) + " " + str(yMax) message = "Maximum extent of input rasters is " + str(gp.extent) + "..."; showPyMessage() #Process: Now build some new grids that have the max extent maxString = "" tempGridList = [] i = 0 tempGridDict = {} for grid in gridsToProcessList: i = i + 1 message = "Normalizing extent of " + grid + "..."; showPyMessage() dsc = gp.describe(grid) if dsc.format != "GRID": #if it aint a grid, make it one tempConvertGrid = gp.createscratchname("","","RASTER",processingWorkspace) gp.extent = "" gp.CopyRaster_management(grid, tempConvertGrid, "", "", "", "", "", "") tempGridList.append(tempConvertGrid) somaExp = "con(isnull(" + tempConvertGrid + "), " + str(tempNoDataValue) + ", " + tempConvertGrid + ")" else: somaExp = "con(isnull(" + grid + "), " + str(tempNoDataValue) + ", " + grid + ")" tempGrid = gp.createscratchname("","","RASTER", processingWorkspace) tempGridList.append(tempGrid) maxString = maxString + tempGrid + "," gp.extent = str(xMin) + " " + str(yMin) + " " + str(xMax) + " " + str(yMax) gp.SingleOutputMapAlgebra_sa(somaExp, tempGrid) #Process: Max, then put back to null message = "Mosaicing - This will take a while..."; showPyMessage() tempGrid = gp.createscratchname("","","RASTER",processingWorkspace) tempGridList.append(tempGrid) if mosaicMethod == "MIN": somaExp = "min(" + maxString + ")" if mosaicMethod == "MAX": somaExp = "max(" + maxString + ")" if mosaicMethod == "MEAN": somaExp = "mean(" + maxString + ")" gp.SingleOutputMapAlgebra_sa(somaExp, tempGrid) somaExp = somaExp = "con(" + tempGrid + " == " + str(tempNoDataValue) + ", setnull(1), " + tempGrid + ")" if outputGrid.split(".")[-1] in (".bmp",".gif",".img",".jpg",".jp2",".png",".tif") or\ ".gdb" in outputGrid or ".mdb" in outputGrid or ".sde" in outputGrid: tempMosaicGrid = gp.createscratchname("","","RASTER", processingWorkspace) tempGridList.append(tempMosaicGrid) gp.SingleOutputMapAlgebra_sa(somaExp, tempMosaicGrid) gp.CopyRaster_management(tempMosaicGrid, outputGrid, "", "", "", "", "", "") else: gp.SingleOutputMapAlgebra_sa(somaExp, outputGrid) #Process: Delete all the temp grids message = "Removing temporay processing grids..."; showPyMessage() for tempGrid in tempGridList: try: gp.Delete_management(tempGrid, "") except: message = "WARNING: Failed to delete processing grid " + str(tempGrid) + "..."; showPyWarning() message = "All done!"; showPyMessage() except: message = "\n*** PYTHON ERRORS *** "; showPyMessage() message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage() message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
Les membres connectés peuvent publier, suivre les mises à jour, et plus encore. Nouveau ici ? Inscrivez-vous gratuitement.
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.