measurementProfile.sml

#
# measurementProfile.sml
# Dan Glasser MI (02Dec03)
# 
# This script takes file input in the format of given by the Geotoolbox's
# Measurement recording feature.  From this file it uses the 'measurement ID'
# (a line number) and Start and End point coordinates.  The two points define
# a line along which DEM sample values will be taken. The sample rate is
# defined according to the Lin/Col Scale  of the raster.
#
# The input file is assumed to have the form:
# (blank line)
# Measurement ID       Length           Angle           Heading          Start X          Start Y           End X            End Y          Centroid X       Centroid Y   
# --------------- ---------------- ---------------- ---------------- ---------------- ---------------- ---------------- ---------------- ---------------- ----------------
# 1                  5643.79691425       8.68975029      81.31024971  635752.11426915 4723968.30723716  641328.90960302 4724835.36975710  638540.51193608 4724401.83849713
# 2                  3617.68325918     -18.43494882     108.43494882  635773.46308871 4728741.99845607  639208.26189042 4727606.30952176  637490.86248957 4728174.15398891
# .
# .
# .
#
# The first 3 lines (the header lines) are ignored.  The other lines are assumed
# to have the values as above separated by a comma or space.  It should be 
# convenient to create this file from measurements taken with the geotoolbox, 
# or user-known coordinates can be used as well.
#
# The results including the Cross-Section Name, Easting, Northing,
# Cumulative Distance, and DEM Z-Value are written to a user-selected file.
# Additional results (for other lines) are concatenated to the end of the file.
#
# declare variables for input file and raster
class RASTER dem;
class FILE inputFile, outputFile;
string defaultPath$ = _context.ScriptDir;
# Get the DEM and the file as input
GetInputRaster(dem);
inputFile = GetInputTextFile("lines.txt", "Please select the profile line file", ".txt");
outputFile = GetOutputTextFile(defaultPath$ + "/results.txt", "Create a file for output", ".txt");
# declare some stuff that we'll need
class Georef georef;
georef = GetLastUsedGeorefObject(dem);
class POINT2D coords;
# declare variables to hold file input
string line$ = "";
numeric x1, y1, x2, y2;
string sectionName$;
# parse the line read from the file to get coord info and section name
proc parseLine(string s$)
{
	sectionName$ = GetToken(line$, ", ", 1, 1);
	x1 = StrToNum( GetToken(line$, ", ", 5, 1) );
	y1 = StrToNum( GetToken(line$, ", ", 6, 1) );
	x2 = StrToNum( GetToken(line$, ", ", 7, 1) );
	y2 = StrToNum( GetToken(line$, ", ", 8, 1) );
}
# given map coordinates, get the DEM value at that point
func getElevationAtPoint(numeric x, numeric y)
{
	local numeric lin, col;
	coords = MapToObject(georef, x, y, dem);
	
	# get the raster line/column number
	col = coords.x;
	lin = coords.y;
	
	return dem[lin, col];
}
# get the header lines and ignore, comment out if header lines don't exist
line$ = fgetline$(inputFile);
line$ = fgetline$(inputFile);
line$ = fgetline$(inputFile);
# get a DEM value every 30m
numeric sampleRate = LinScale(dem), samples = 0;
numeric xDist = 0, yDist = 0, m = 0, b = 0, deltaX;
numeric x, y;
# loop through the rest of the lines to get coordinate info
while (feof(inputFile)==0)
{
	sampleRate = abs(sampleRate); # re-positivify the sample rate
	# read and handle a line at a time
	line$ = fgetline$(inputFile);
	parseLine(line$);
	# print the header line for the output file
	fprint(outputFile, "Cross-Section,", "Easting,", "Northing,", "Distance,", "Z-Value");
	
	# get the absolute distance from point A to point B
	numeric distance = 0;
	xDist = abs(x2 - x1);
	yDist = abs(y2 - y1);
	samples = round( sqrt( sqr(xDist) + sqr(yDist) ) / sampleRate );
	
	# handle the non-vertical line case
	if (x1 != x2)
	{
		# compute the slope (m) of the line ( y = m x + b ) 	
		m = (y2 - y1) / (x2 - x1);
		b = y1 - m * x1;
		x = x1;
		deltaX = (x2 - x1) / samples;	# x stepsize		
		
		numeric j;
		for j = 0 to samples
		{
			y = m * x + b;		
			fprint(outputFile, sectionName$, x, y, distance, getElevationAtPoint(x, y));
			x = x + deltaX;
			distance = distance + sampleRate;
		}
	}
	# handle the vertical line case
	else
	{		
		if(y2<y1)
		{
			sampleRate = -sampleRate;
		}
		x = x1;
		y = y1;
	
		numeric j;
		for j = 0 to samples
		{
			fprint(outputFile, sectionName$, x, y, distance, getElevationAtPoint(x, y));
			y = y + sampleRate;
			distance = distance + sampleRate;
		}	
	}
}
# close handles
CloseRaster(dem);
fclose(inputFile);
fclose(outputFile);