coastal.sml

########################################################
#
#	Transect and Rates SML script
#
#
#	Written 8 Mar 1998
#	David Williss
#	MicroImages, Inc
#
#	Based on a pair of DSAS programs.
#
########################################################
#
#	Limits --
#		A few limits which are hard-coded.  
#		Increase as necessary
#
numeric maxTrans = 10000;		# Maximum number of transects
numeric maxShores = 20;			# Maximum number of shoreline vectors

########################################################
#
#	Defaults --
#		These are initial default values.  
#		They will be read out of the tntproc.ini file.  The
#		last parameter to IniReadNumber is the default value
#		to use if not found in the ini file
#
numeric transInterval, maxdist;
string outfilebase$;
transInterval = IniReadNumber("transect", "interval", 100);
maxdist = IniReadNumber("transect", "maxdist", 250);
outfilebase$ = FileNameGetPath(_context.Filename);
outfilebase$ = IniReadString("transect", "outfilebase", outfilebase$);

########################################################
#
#	Global Variables
#
array numeric  transx[maxTrans];
array numeric transy[maxTrans];
array numeric transa[maxTrans];
numeric numtrans, numshores;
array numeric shorex[maxShores], shorey[maxShores];
array numeric vectdate[maxShores], juldate[maxShores], years[maxShores];
array numeric distToBase[maxShores], distFromShore1[maxShores], years[maxShores];
array numeric basex[10], basey[10];	# Automatically resized anyway
class FILE locfile, julfile, distfile, ratefile;
class GRE_GROUP group;
class GRE_LAYER_VECTOR voverlay, layer;
class GRE_LAYER translayer;
class GRE_VIEW view;
class GRE_LAYER_RASTER bglayer;
class XmForm form;
class PUSHBUTTONITEM okItem, cancelItem;
class XmForm button_row;
class XmSeparator sep;
class PromptNum pinterval, pmaxdist;
class STATUSCONTEXT status;
class DATABASE database;
class DBTABLEINFO ratetable;
class DBFIELDINFO transfield, field;
class GRE_LAYERMANAGER controls;
class XmPushButton transButton;
vector Transects;
raster BackgroundRaster;
vector BaseLineVector;
string prompt$, fname$, name$, lastdate$;
numeric inode, year, month, date;
numeric line, numbasepts, flipbaseline;
numeric x1, y1, x2, y2;
numeric hits, misses, numint, sx, sy, rates;
########################################################

proc DoTransSeg(x0,y0, x1,y1) 
{
	local numeric theta, theta90, r, dx, dy, x, y, dist;

	dx = x1 - x0;
	dy = y1 - y0;
	theta = atan2(dy, dx);	# Angle of this segment
	theta90 = theta + pi/2;			# angle perpendicular to segment
	dist = sqrt(dx*dx + dy*dy);
	for r = 0 to dist step transInterval {
		numtrans += 1;
		transx[numtrans] = x0 + r*cos(theta);
		transy[numtrans] = y0 + r*sin(theta);
		transa[numtrans] = theta90;
		}
	if (r != dist) {
		numtrans += 1;
		transx[numtrans] = x1;
		transy[numtrans] = y1;
		transa[numtrans] = theta90;
		}
	return;
   }

########################################################
#
#	Do the corner between two segments. 
#		Segment A is from (x0,y0) to (x1,y1)
#		Segment B is from (x1,y1) to (x2,y2)
#
#	The transect from the end of segment A is already in the
#	list.  The transect from the beginning of B will be added
#	by DoTransSeg
#	Compute a transect at point (x1,y1) which bisects the
#	angle betweent the two transects.
#
proc DoCorner(x0,y0, x1,y1, x2,y2) {
	local numeric dxA, dxB, dyA, dyB, thetaA, thetaB, theta90A, theta90B, theta90AB;

	dxA = x1 - x0;
	dyA = y1 - y0;
	dxB = x2 - x1;
	dyB = y2 - y1;
	thetaA = atan2(dyA, dxA);	# Angle of segment A
	thetaB = atan2(dyB, dxB);	# Angle of segment B
	theta90A = thetaA + pi/2;
	theta90B = thetaB + pi/2;
	# Angles greater than 2*pi screw up the bisection
	if (theta90A > 2 * pi) { 
		theta90A -= 2 * pi; 
		}
	if (theta90B > 2 * pi) { 
		theta90B -= 2 * pi; 
		}
	if (abs(theta90A - theta90B) > pi) {
		theta90A -= 2 * pi;	# Doesn't really matter which one we subtract from
		}
	theta90AB = (theta90A + theta90B) / 2;
	numtrans += 1;

	transx[numtrans] = x1;
	transy[numtrans] = y1;
	transa[numtrans] = theta90AB;
	return;
	}

########################################################
#
#		Compute the endpoint rate dr/dt.
#		The dates are assumed to be julian dates
#		return value is in meters/year.
#
func EndPointRate(dist0, dist1, date0, date1) {
	return ((dist1 - dist0) / ((date1 - date0) / 365.25));
	}


#	Called when user closes layer controls.  
#	Close the view window.  This will cause
#	the popdown callback to be called
proc cbQuit(class widget widget) {
	DialogClose(form);
	}

#	Called when the dialog actually pops down.  Free
#	and destroy everything.
proc cbPopdown(class widget widget) {
	ViewDestroy(view);
	GroupDestroy(group);
	DestroyWidget(form);
	Exit();
	}


#######################################################
#								MAIN 
#######################################################

#
#	Lets assume that the baseline vector contains lines whos
#	verticies are the segment points.
#
if (_context.ScriptIsRVCObject) {
	# Script is being run from an RVC file which has all the data it needs
	# This is for demo purposes
	OpenVector(BaseLineVector, _context.Filename, "BaseLine");
	numshores = OpenInputVectorList(ShoreVectorList, 
				_context.Filename, "_1844",
				_context.Filename, "_1942",
				_context.Filename, "_1976"
				);
	}
else {
	GetInputVector(BaseLineVector);
	numshores = GetInputVectorList(ShoreVectorList, "Select shoreline vectors in cronological order");
	}

class XmForm pform;
numeric cancel;

proc cbPromptsOK() {
	DialogClose(pform);
	cancel = 0;
	}

proc cbPromptsCancel() {
	DialogClose(pform);
	cancel = 1;
	}

proc PromptForValues() {
	pform = CreateModalFormDialog("Parameters");
	okItem = CreatePushButtonItem("OK", cbPromptsOK);
	cancelItem = CreatePushButtonItem("Cancel", cbPromptsCancel);
	button_row = CreateButtonRow(pform, okItem, cancelItem);
	button_row.BottomWidget = pform;

	sep = CreateHorizontalSeparator(pform);
	sep.LeftWidget = pform;
	sep.RightWidget = pform;
	sep.BottomWidget = button_row;

	pinterval = CreatePromptNum(pform, "Interval:", 5, 0, transInterval, 0);
	pinterval.BottomWidget = sep;

	pmaxdist = CreatePromptNum(pform, "Max Dist:", 5, 0, maxdist, 0);
	pmaxdist.BottomWidget = pinterval;
	AlignWidgets(pinterval.label, pmaxdist.label);

	DialogOpen(pform);
	DialogWaitForClose(pform);
	if (!cancel) {
		transInterval = pinterval.value;
		maxdist = pmaxdist.value;
		# Write them back out to the ini file for next time
		IniWriteNumber("transect", "interval", transInterval);
		IniWriteNumber("transect", "maxdist", maxdist);
		}
	DestroyWidget(pform);
	}

proc cbComputeTransects() {
	local i;

	PromptForValues();
	if (cancel) {
		return;
		}
	#------------------------------------------------------------
	#	Prompt for output file location
	#		I'd like to make this use the standard get file
	#		dialog but I don't have that callable from within
	#		SML yet.
	#
	prompt$ = "Select output .loc file.  Other files will\n"
			+ "be created in the same directory by changing\n"
			+ "the file extention";
	GetOutputVector(Transects, "Network,VectorToolkit");
	fname$ = Transects.$info.Filename;
	#fname$ = GetOutputFileName(outfilebase$ + ".loc", prompt$, ".loc");
	outfilebase$ = FileNameGetPath(fname$) + FileNameGetName(fname$);
	IniWriteString("transect", "outfilebase", outfilebase$);	# remember it for next time

	#------------------------------------------------------------
	#	Prompt operator for the dates of the shoreline vectors.  
	#	Genreate the default by assuming the vector names are of
	#	the form _YYYY_MM
	#
	numeric j;
	for j = 1 to numshores {
		fname$ = GetObjectFileName(ShoreVectorList[j]);
		inode = GetObjectNumber(ShoreVectorList[j]);
		name$ = GetObjectName(fname$, inode);
		year = StrToNum(mid$(name$, 2, 4));
		month = StrToNum(mid$(name$, 7, 2));
		if (month < 1) month = 1;
		date = year * 10000 + month * 100 + 1;
		if (! _context.ScriptIsRVCObject) {
			# Don't bother prompting if this is only a demo.  Just assume Jan 1
			prompt$ = "Enter date for vector " + name$ + " as YYYYMMDD";
			date = PopupNum(prompt$, date);
			}
		vectdate[j] = DateToJulian(date);
		}


	status = StatusContextCreate(view.StatusHandle);
	#------------------------------------------------------------
	#
	#	Initialize vector toolkit routines for intersection
	#	finding on the shoreline vectors
	#
	for j = 1 to numshores {
		status.Message = sprintf("Preparing vector %d for finding intersections\n", j);
		VectorToolkitInit(ShoreVectorList[j]);
		}

	#------------------------------------------------------------
	#
	#	Compute the transects
	#
	numtrans = 0;
	clear();
	status.Message =  "Computing transects...";
	for line = 1 to NumVectorLines(BaseLineVector) {
		numbasepts = GetVectorLinePointList(BaseLineVector, basex, basey, line);
		if (flipbaseline) {
			for i = 1 to numbasepts - 1 {
				if (i > 1) {
					DoCorner(basex[i+1], basey[i+1], basex[i], basey[i], basex[i-1], basey[i-1]);
					}
				DoTransSeg(basex[i+1], basey[i+1], basex[i], basey[i]);
				}
			}
		else {
			for i = 1 to numbasepts - 1 {
				if (i > 1) {
					DoCorner(basex[i-1], basey[i-1], basex[i], basey[i], basex[i+1], basey[i+1]);
					}
				DoTransSeg(basex[i], basey[i], basex[i+1], basey[i+1]);
				}
			}
		}

	#------------------------------------------------------------
	#
	#	Create a vector object with the trancects so we can overlay them
	#	Would like to make this use network topology.  That way, transects
	#	which cross each other won't be broken up into multiple lines.  This
	#	insures that transect 1000 is line 1000 and will make it easier to
	#	attach database records later
	#
	status.Message = "Creating transect vector...";
	#DeleteFile(outfilebase$ + ".rvc");
	#CreateVector(Transects, outfilebase$ + ".rvc", "Transects", "Generated by SML", "Network");
	for i = 1 to numtrans {
		StatusSetBar(status, i, numtrans);
		x1 = transx[i];
		y1 = transy[i];
		x2 = x1 + maxdist * cos(transa[i]);
		y2 = y1 + maxdist * sin(transa[i]);
		VectorAddTwoPointLine(Transects, x1, y1, x2, y2);
		}
	StatusSetBar(status);
	status.Message = "Updating transect vector topology...";
	CopySubobjects(BaseLineVector, Transects, "georef");
	VectorValidate(Transects);
	CloseVector(Transects);

	database = OpenVectorLineDatabase(Transects);
	ratetable = TableCreate(database, "Rates", "Erosion rates");
	ratetable.ImpliedAttachment = 1;
	transfield = TableAddFieldInteger(Transects.line.Rates, "Transect", 5);
	field = TableAddFieldFloat(Transects.line.Rates, "EPR", 8, 2);
	field.UnitType = "velocity";
	field.Units = "m/yr";
	field = TableAddFieldFloat(Transects.line.Rates, "Mean", 8, 2);
	field.UnitType = "velocity";
	field.Units = "m/yr";
	TableAddFieldFloat(Transects.line.Rates, "StdDev", 8, 2);
	TableAddFieldFloat(Transects.line.Rates, "Variance", 8, 2);
	TableAddFieldFloat(Transects.line.Rates, "Regression", 8, 2);
	TableAddFieldFloat(Transects.line.Rates, "Jackknife", 8, 2);

	if (! _context.ScriptIsRVCObject) {
		# Don't create text files for the demo version
		# If we don't open the files then fprintf or close,
		# SML is kind enough to just ignore us (C programs
		# die if we did that)
		julfile = fopen(outfilebase$ + ".jul", "w");
		distfile = fopen(outfilebase$ + ".dis", "w");
		ratefile = fopen(outfilebase$ + ".rat", "w");
		locfile = fopen(outfilebase$ + ".loc", "w");
		}

	status.Message = "Creating .loc file...";

	#------------------------------------------------------------
	#
	#	Output the .loc file
	#
	for i = 1 to numtrans {
		StatusSetBar(status, i, numtrans);
		fprintf(locfile, "%f,%f,%d\n", transx[i], transy[i], i);
		}


	#------------------------------------------------------------
	#	This is the biggie!  Find the intersections and compute
	#	the rates.
	#
	status.Message = "Finding shoreline-transect intersections";

	hits = 0;
	misses = 0;
	for i = 1 to numtrans {
		StatusSetBar(status, i, numtrans);

		#------------------------------------------------------------
		# Find intersection of vector transx[i], transy[i], transa[i]
		# with each line.  
		#
		#	Write the intersections to the .jul file.
		#
		#	Point to ask Randy: do you want this as one line per transect like
		#	the original or one line per intersection as requested for the .dis
		#	file?
		#
		numint = 0;
		for j = 1 to numshores {
			if (VectorLineRayIntersection(ShoreVectorList[j], transx[i], transy[i], transa[i], maxdist, sx, sy)) {
				fprintf(julfile, "%4d %10.2f %10.2f %10.2f %10.2f %7d\n", i, transx[i], transy[i], sx, sy, juldate[j]);
				numint += 1;
				shorex[numint] = sx;
				shorey[numint] = sy;
				distToBase[numint] = sqrt((transx[i] - sx) ^ 2 + (transy[i] - sy) ^ 2);
				juldate[numint] = vectdate[j];
				}
			}
		if (numint) {
			hits += 1;
			}
		else {
			misses += 1;
			}

		#------------------------------------------------------------
		#	Output the line in the .dist file
		#
		#	For each intersection, also compute the distance from this
		#	intersection and the first intersection, and also the
		#	difference in time (in years).  These will obviously be
		#	0 for the first intersection.
		#
		string thisdate$, lastdate;
		numeric lastdist = 0, delta;
		for j = 1 to numint {
			thisdate$ = DateToString(JulianToDate(juldate[j]), "%m/%Y");
			if (j > 1) {
				delta = distToBase[j] - lastdist;
				fprintf(distfile, "%d,%s,%s,%f\n", i, lastdate$, thisdate$, delta);
				}
			distFromShore1[j] = sqrt((shorex[j] - shorex[1])^2 + (shorey[j] - shorey[1])^2);
			if (distToBase[1] > distToBase[j]) distFromShore1[j] = -distFromShore1[j];
			years[j] = (juldate[j] - juldate[1]) / 365.25;

			lastdist = distToBase[j];
			lastdate$ = thisdate$;
			}

		#------------------------------------------------------------
		#
		#	Compute endpoint rate
		#
		fprintf(ratefile, "%d,", i);	# transect number
		Transects.line[i].Rates.Transect = i;
		if (numint > 1) {
			numeric epr = EndPointRate(distToBase[numint], distToBase[1], juldate[numint], juldate[1]);
			fprintf(ratefile, "%.2f", epr);
			Transects.line[i].Rates.EPR = epr;
			}

		#------------------------------------------------------------
		#
		#	Compute average of rates.  Need at least 3 points to do this
		#		This uses a feature of SML that's been there since the beginning
		#		but I don't know if it's documented.  The rates variable is set to
		#		the empty set ( { } ).  Then each endpoint rate is computed and
		#		the is added to the set.  When the loops are done, we have a set
		#		of all the rates, which can be passed to SetMean, SetSD, SetVariance
		#		or any of the other Set functions.
		#
		numeric k;
		if (numint >= 2) {
			rates = {};
			for j = 1 to numint - 1 {
				for k = j+1 to numint {
					epr = EndPointRate(distToBase[j], distToBase[k], juldate[j], juldate[k]);
					rates = rates + { epr };
					}
				}
			Transects.line[i].Rates.Mean = SetMean(rates);
			fprintf(ratefile, ",%.2f", SetMean(rates));
			if (numint > 2) {
				fprintf(ratefile, ",%.2f,%.2f", SetSD(rates), SetVariance(rates));
				Transects.line[i].Rates.StdDev = SetSD(rates);
				Transects.line[i].Rates.Variance = SetVariance(rates);
				}
			else {
				fprintf(ratefile, ",,");
				}
			}
		else {
			fprintf(ratefile, ",,,");
			}

		#------------------------------------------------------------
		#
		#	Standard linear regression.  Compute slope where x is years and
		#	y is distance in meters.
		#
		numeric slope, intercept;
		if (numint > 1) {
			LinearRegression(years, distFromShore1, numint, slope, intercept);
			fprintf(ratefile, ",%.2f", slope);
			Transects.line[i].Rates.Regression = slope;
			}
		else {
			fprintf(ratefile, ",");
			}

		#------------------------------------------------------------
		#
		#	Compute jackknife regression.   The jackknife regression is the 
		#		mean of the regression run numint times, each time omitting 
		#		one of the points.  The LinearRegression function has an
		#		optional parameter that makes this easier.  The parameter 
		#		after the intercept tells it which point to ignore. 
		#
		numeric sum;
		if (numint > 2) {
			sum = 0;
			for j = 1 to numint {
				LinearRegression(years, distFromShore1, numint, slope, intercept, j);
				sum += slope;
				}
			fprintf(ratefile, ",%.2f", sum/numint);
			Transects.line[i].Rates.Jackknife = sum/numint;
			}
		else {
			fprintf(ratefile, ",");
			}

		fprintf(ratefile, "\n");
		}

	StatusSetBar(status);	# Clear it
	status.Message = "Cleaning up...";

	for i = 1 to  numshores {
		CloseVector(ShoreVectorList[i]);
		}

	fclose(julfile);
	fclose(distfile);
	fclose(ratefile);
	fclose(locfile);
	StatusContextDestroy(status);

	if (voverlay) {
		LayerDestroy(voverlay);
		}
	voverlay = GroupQuickAddVectorVar(group, Transects);
	voverlay.Line.NormalStyle.Color.name = "blue";
	voverlay.Line.DataTip.Field = FieldGetInfoByName(Transects.line.Rates, "Mean");
	ViewRedraw(view);
	}

#############################################################
############################ MAIN ###########################
#############################################################

	numeric height = .8 * _context.DisplayInfo.Height;
	numeric width = .75 * height;
	form = CreateFormDialog("Transect intersections");
	WidgetAddCallback(form.Shell.PopdownCallback, cbPopdown);

	group = DispCreate2DGroup();
	if (_context.ScriptIsRVCObject) {
		if (RasterExists(_context.Filename, "KENTIS")) {
			OpenRaster(BackgroundRaster, _context.Filename, "KENTIS");
			bglayer = GroupQuickAddRasterVar(group, BackgroundRaster);
			bglayer.ShowDataTips = 0;
			width = height * (BackgroundRaster.$info.NumCols /
								   BackgroundRaster.$info.NumLins);
			}
		}
	numeric i;
	for i = 1 to numshores {
		layer = GroupQuickAddVectorVar(group, ShoreVectorList[i]);
		}
	layer = GroupQuickAddVectorVar(group, BaseLineVector);
	layer.Line.NormalStyle.Color.name = "blue";

	view = GroupCreateView(group, form, "view", height, width);
	view.form.TopWidget = form;
	view.form.LeftWidget = form;
	view.form.RightWidget = form;
	view.form.BottomWidget = form;
	view.BackgroundColor.Red = 80;
	view.BackgroundColor.Green = 80;
	view.BackgroundColor.Blue = 80;
	ViewCreateSelectTool(view);
	ViewCreateMeasureTool(view);
	ViewAddStandardTools(view);
	ViewAddToolIcons(view);
	transButton = CreatePushButton(view.IconBar, "Create Transects");
	WidgetAddCallback(transButton.ActivateCallback, cbComputeTransects);
	DialogOpen(form);
	controls = GroupOpenLayerManagerWindow(form, group);
	DispAddCallback(controls.CloseCallback, cbQuit);
	ViewRedraw(view);
	WaitForExit();