Dicing large mosaics

Question

What is the best way to take a large mosaiced raster dataset (e.g. from http://viewer.nationalmap.gov/viewer/) and tile it to USGS 3.75 minute quarter quad specs? The problem I am having is creating 300m overlapping borders and applying USGS naming conventions to the output. Thanks in advance.

Answer 1

What software do you have? A real quick solution is Global Mapper (which costs $350+). If you have this it is simply file-->export raster --> choose file type and then in the 3rd tab (gridding) you can import a vector file with your quad or define it. It is also possible in arcpy and Grass (will add details if you want but it all depends on what software you want to utilize)

=== Here is what I use in python / arcpy.

You will need to modify the logic of the naming convention.

# Clips raster to polygon
# Outputs csv with original and new names
# Naming convention specific to current fileset and MUST be re-written
# ovr is the ESRI Raster Pyramid

#Licence: Creative Commons
#Created by: George Corea; georgec@atgis.com.au, coreagc@gmail.com


import arcpy, glob, os, sys, arcgisscripting
from arcpy import env, mapping

# SET UP YOUR WORKING ENVIRONMENT, INPUT and OUTPUT

path = os.getcwd()
env.workspace = path
os.chdir(path)
env.overwriteOutput = True


RasterInputType = ['tif'] # Type in 'tif, ecw' etc the type of Rasters which need to be processed
RasterOutputType = 'tif' #BIL, BIP, BMP, BSQ, DAT, GIF, GRID, IMG, JPG, JP2, PNG, TIFF, or any geodatabase raster dataset. ECW is not supported as output
rasterType = "INTEGER"
spatial_reference = arcpy.SpatialReference("C:\\Python26\\GDA_1994_MGA_Zone_55.prj") # point to the required projection
AOI = 'AOI_edited_region.shp' # type in the name or path to and name of the clipping extent

#END COMMON set up parameters except for file name.

print 'Reading files from ' + path

try:
    os.mkdir (RasterOutputType+"_clip")

except:
    os.remove (RasterOutputType+'\\*')

x=0


for Raster in RasterInputType:
    FileList = glob.glob("*." + Raster)
    FileListCount = len(FileList)
    print "Processing "+str(FileListCount)+" files in directory."
    for File in FileList:
        # The folowing describes the input and output filenames. \n
        # This will need to be modified for datasets other than the current lidar
        PositionDot =File.find('.')
        Position1 = File.find('000_')
        Position2 = File.find('000.')
        Position1b = Position1-4
        Position2b = Position2-4
        #ext =File[PositionDot:]
        Name = File[0:PositionDot]
        Name1 = File[0:Position1] # Character string to retain from xth character to yth character
        Name2 = File[Position2b:Position2]
        FullName = File[0:]   
        inFile = path+"\\"+File
        outRasterName = Name1+"-"+Name2+"_"+str(x)+'.'+RasterOutputType
        outRaster = path+'\\'+RasterOutputType+'_clip\\'+outRasterName
        print inFile, outRasterName, outRaster
        # END describing in/out files

        # Execute Conversion function
        arcpy.Clip_management(inFile, "#", outRaster, AOI, "0", "ClippingGeometry")
        # Output file
        f = open('RasterClip_OUTPUT.txt', 'a')
        f.write(str(Name)+","+str(outRasterName)+"\n")
        f.close()
        print "Using "+inFile+" CREATED "+ outRaster
        outfile= path+'\\'+RasterOutputType+'_clip\\'+FullName
        #File.close
        #os.rename(File, outfile)
        x=x+1

print "Processed and created "+ str(x)+ "files and RasterClip_OUTPUT.txt"