# 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()
Angemeldete Mitglieder können Beiträge verfassen, Updates folgen und mehr. Neu hier? Registriere ein kostenloses Konto.
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.