|
Although this Federal Geographic Data Committee-compliant metadata file is intended to document the data set in nonproprietary form, as well as in ArcInfo format, this metadata file may include some ArcInfo-specific terminology.
Ran Python script, "Bathy Processor v1.5", designed to filter the bathymetry text file
****Bathy_Proc (Python Script)****
#!/usr/bin/python
# Process bathymetry for P1-06-MB
# Once this stabalizes, I will code into either a batch file or something more robust
import os, glob, sys, time
# Configure Directories
SCRIPT_DIR = "E:/subx_filters"
TXT_DIR = "D:/S205SC/Bathy/TXT"
GRID_DIR = "D:/S205SC/Bathy/asciigrids"
# Trackline Filtering Parameters
MAX_SWATH_MULTIPLE = 4.0 # Trim depth samples > MAX * Water Depth (trims far-field noise)
MIN_SWATH_MULTIPLE = 0.0 # Trim depth samples < MIN * Water Depth (trims nadir noise)
ROLL_BEAMS = 150 # Number of beams to average to eliminate roll artifacts
NOISE_BEAMS = 300 # Number of beams to average to eliminate along-track noise
# Output Grid Parameters
GRID_RESOLUTION = 1 # Output grid resolution
SEARCH_RADIUS = GRID_RESOLUTION + 5 # Search radius for samples while gridding
# Get the files for processing
files = glob.glob("%s/*txt" % TXT_DIR)
# Get the total size of the input files to estimate processing progress
bytes_total = sum([os.path.getsize(x) for x in files])
bytes_sofar = 0
for i, file in enumerate(files):
time_start = time.time()
(indir, infile) = os.path.split(file)
(basename, ext) = os.path.splitext(infile)
print
print "-------------------------------------------------------------------------"
print "Begin Processing line: %s (%d of %d)" % (basename, i+1, len(files))
print
if not os.path.exists(basename[0:9]):
print "Step 1 of 8: converting file to binary format..."
outfile = basename + ".bin"
os.system('python %s/seaconvert.py -bnxzad -o %s %s' % (SCRIPT_DIR, outfile, file))
print "Step 2 of 8: trimming beams..."
infile = outfile
outfile = basename + ".bin.flt"
os.system('python -u %s/trim.py -d%f -n%f %s > %s' % (SCRIPT_DIR, MAX_SWATH_MULTIPLE, MIN_SWATH_MULTIPLE, infile, outfile))
os.remove(infile)
print "Step 3 of 8: removing roll artifacts..."
infile = outfile
outfile = basename + ".bin.flt.roll"
os.system('python -u %s/deroll.py %s %d > %s' % (SCRIPT_DIR, infile, ROLL_BEAMS, outfile))
os.remove(infile)
print "Step 4 of 8: removing system noise..."
infile = outfile
outfile = basename + ".bin.flt.roll.noise"
os.system('python -u %s/denoise.py %s %d > %s' % (SCRIPT_DIR, infile, NOISE_BEAMS, outfile))
os.remove(infile)
print "Step 5 of 8: calculate spatial extent of line..."
infile = outfile
outfile = basename + ".region"
os.system('gmtset D_FORMAT %.2f')
os.system('gmtconvert %s -bi6 -F1,2 -bo > temp1' % infile)
os.system('minmax -bi2 -I%d temp1 > %s' % (GRID_RESOLUTION, outfile))
os.remove('temp1')
print "Step 6 of 8: gridding line..."
Rfile = basename + ".region"
infile = basename + ".bin.flt.roll.noise"
outfile = basename + ".nc"
region = open(Rfile).readline().strip()
os.system('gmtconvert %s -bi6 -F1,2,3 -bo > temp1' % (infile))
os.system('blockmean -I%d %s -C -bi3 -bo temp1 > temp2' % (GRID_RESOLUTION, region))
os.remove('temp1')
os.system('nearneighbor -G%s -I%d -bi3 -N4 -R -S%d temp2' % (outfile, GRID_RESOLUTION, SEARCH_RADIUS))
os.remove('temp2')
os.remove(infile)
os.remove(Rfile)
print "Step 7 of 8: convert line to ESRI ASCII Grid format..."
infile = basename + ".nc"
outfile = GRID_DIR + "/" + basename + ".asc"
os.system('grdmath %s 100 MUL RINT = temp.nc' % (infile))
os.system('grd2xyz temp.nc -E > %s' % (outfile))
os.remove('temp.nc')
os.remove(infile)
print "Step 8 of 8: convert line to ESRI Raster format..."
infile = GRID_DIR + "/" + basename + ".asc"
outfile = GRID_DIR + "/" + basename
#os.system('ascii2arc.py %s %s' % (infile, outfile))
os.system("ascii_to_arcraster_project_pyramids.py %s %s" % (infile, outfile))
os.remove(infile)
# Calculate Estimated time remaining
time_stop = time.time()
if time_stop == time_start:
time_stop = time_start + 1
bytes_file = os.path.getsize(file)
bytes_sofar = bytes_sofar + bytes_file
bytes_permin = float(bytes_file) / (time_stop - time_start) * 60.0
bytes_remaining = bytes_total - bytes_sofar
time_remaining = bytes_remaining / bytes_permin
percent_complete = float(bytes_sofar) / float(bytes_total) * 100
# End this line
print
print "Line Completed. Completed: %.0f%%, Est. Time Remaining: %.0f minutes" % (percent_complete, time_remaining)
print "-------------------------------------------------------------------------"
print
Ran Python script, "ascii2arc", designed to convert ASCII grids to ESRI ArcGrid fromat
****ascii2arc Script)****
# ---------------------------------------------------------------------------
# ascii_to_arcraster_project_pyramids.py
# Created on: Mon Oct 02 2006 01:45:01 PM
# Nadine Golden
# Usage: ascii2arc.py <temp1_txt> <Output_raster>
# ---------------------------------------------------------------------------
# Import system modules
import glob, sys, string, os, win32com.client
# Create the Geoprocessor object
gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")
# Load required toolboxes...
gp.SetProduct("ArcInfo")
gp.AddToolbox("E:/Mapping/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
gp.AddToolbox("E:/Mapping/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
# Get script arguments...
if len(sys.argv) < 3:
print >> sys.stderr, "ascii_to_arcraster_project_pyramids.py: error1: not enough arguments"
sys.exit(1)
# Convert ascii grid to arc grid
ingrid = sys.argv[1]
outgrid = sys.argv[2]
(tail, head) = os.path.split(ingrid)
(root, ext) = os.path.splitext(head)
if len(root) > 13:
print >> sys.stderr, "ascii2arc.py: error 1: Grid name exceeds 13 characters: %s" % root
sys.exit(1)
gp.Workspace = tail
gp.ScratchWorkspace = tail
#print "Workspace:", gp.Workspace
tmpgrid = tail + "/temp"
# Check to see if temp already exists and delete it if necessary
if os.path.exists(tmpgrid):
gp.Delete(tmpgrid)
#print "ingrid1: %s outgrid1: %s" % (ingrid, tmpgrid)
gp.ASCIIToRaster_conversion(ingrid, tmpgrid, "INTEGER")
# Process: Project Raster...
ingrid = tmpgrid
outgrid = tail + "/" + root
if os.path.exists(outgrid):
gp.Delete(outgrid)
gp.ProjectRaster_management(ingrid, outgrid,
"PROJCS['WGS_1984_UTM_Zone_10N',
GEOGCS['GCS_WGS_1984',
DATUM['D_WGS_1984',
SPHEROID['WGS_1984',6378137.0,298.257223563]],
PRIMEM['Greenwich',0.0],
UNIT['Degree',0.0174532925199433]],
PROJECTION['Transverse_Mercator'],
PARAMETER['False_Easting',500000.0],
PARAMETER['False_Northing',0.0],
PARAMETER['Central_Meridian',-123.0],
PARAMETER['Scale_Factor',0.9996],
PARAMETER['Latitude_Of_Origin',0.0],
UNIT['Meter',1.0]];
-10000 -10000 100000;0 100000;0 100000",
"NEAREST", "1")
# Process: Build Pyramids...
gp.BuildPyramids_management(outgrid)
# Clean up the mess
if os.path.exists(tmpgrid):
gp.Delete(tmpgrid)
Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.