##################################################################################
#
# StrikeDipTool.sml (Tool Script)
#
##################################################################################
#
# Started: 3 November 2003
# Updated: 19 June 2007 to work in overhead stereo mode in the main view.
#
# Authors:
# Dave Breitwisch
# Randy Smith
# Dan Glasser
#
# Version Nov. 2007
# Requires TNTmips 2007:73 or later
# Revisions:
# Produces correct attitude values for DEMs with geographic coordinates by
# using appropriate UTM zone coordinates for attitude computations.
# Automatically applies raster scale value to allow scaled DEMs with values in feet or decimeters.
# Strike azimuth and dip direction are automatically adjusted for the local
# grid angle to north at each point.
#
# Version 24 March 2008
# Fixes a problem with the user-defined function to return angle to north
# that is used to correct the strike azimuth
##################################################################################
##################################################################################
#
# Purpose:
#
# The purpose of this Tool Script is to give a geologist a streamlined method
# for determining the Strike and Dip values from an accurate DEM raster and
# overlaid arial imagery. The accompanying CartoScript displays the added
# points using the appropriate attitude symbols. The script creates two
# vector objects: a StrikeDipVector to store point elements with attached
# attitude values in the associated point database, and an OutcropTrace
# object to store outcrop traces for each point (see below). The script creates
# and displays a dialog to show the computed attitude values and allow
# designation of overturned beds.
#
# Assumptions:
#
# The DEM is the bottom layer in the active group.
#
# Description:
#
# This Tool Script has four different modes that the geologist can use:
#
# ADD - While viewing the reference image, place 3 points (three vertices of a
# triangle) on a bedding plane by left clicking on the screen using the
# provided polygon tool. Adjust the vertices if needed using the
# Line/Polygon Edit Controls. Right-click or press [Apply] on the
# Edit Controls: x = easting and y = northing values are registered for
# each point from the group map coordinates and z = elevation is registered
# for each point from the corresponding cell of the elevation raster that
# is the first layer in the group. If the xyz values pass a validity check,
# these three noncolinear points are used to compute the strike azimuth,
# dip angle, and dip azimuth of the plane. These values are shown in
# the Strike-Dip Toolscript Controls window, which also includes a
# toggle button to indicate if the bed is overturned.
#
# The Tool Script also creates and displays an outcrop trace in the vicinity
# of the attitude measurement. This trace is the computed line of intersection
# between the measured plane and the surface represented by the DEM.
# Verify that the attitude values appear reasonable and that the computed
# outcrop trace parallels the actual outcrop trace in the image. (The
# Line/Polygon tool remains active to allow you to adjust the points and
# reapply to recompute a new attitude measurement and outcrop trace if
# needed.)
#
# When you are satisfied, press the Accept button on the Tool Script
# Controls window. A point is added to the StrikeDip vector object at the
# center of the triangle and the attitude values are written to an attached
# record in the associated point element database. The point is styled
# automatically with the appropriately-oriented and labeled attitude symbol
# using the accompanying CartoScript. The CartoScript must be in the
# same directory as this Tool Script (or, if the Tool Script is saved with
# a Group or Layout, in the same directory as the Project File
# containing that Group or Layout). The outcrop trace line is also saved
# in the separate OutcropTrace vector object along with a specified color.
# The original three points are also saved in this color with the trace line.
#
# DELETE - Select a previously-added Strike/Dip point and press [Accept] to delete it
# along with the outcrop trace pattern that is associated with it.
#
# EDIT - Select a previously-added Strike/Dip point and press [Accept] to
# edit it. The original triangle and outcrop trace are reloaded to
# enable vertex locations to be moved using the Line/Polygon tool as
# in Add mode; Apply or right-click to recompute a new attitude and
# outcrop tract. Press [Accept] on the Tool Script Control window to
# save the edited point and outcrop trace.
#
# VIEW - Select a Strike/Dip point and use the Accept button to toggle whether
# the outcrop trace pattern associated with that point is displayed.
# Both the trace pattern and the three points are displayed in their
# saved color.
#
# Definitions:
#
# Strike Angle = azimuth of the line of intersection of the plane with the horizontal.
# Note: The Strike Angle value conforms to the right hand rule, where the strike
# angle is in the direction along the strike line for which the dip
# direction is to the right.
#
# Dip Angle = maximum vertical angle of the plane measured downward from the horizontal.
#
# Dip Direction = azimuth of the dip vector, by definition is perpendicular to Strike.
#
# Linear equation of the plane = a*x + b*y c*z + d = 0
#
# Note: The coordinates of the three selected points are used to compute the
# coefficient values of a, b and c. d is then solved mathematically.
#
##################################################################################
##################################################################################
#
# Standard ToolScript Comments:
#
# The following symbols are predefined
# class GRE_VIEW View {use to access the view the tool script is attached to}
# class GRE_GROUP Group {use to access the group being viewed if the script is run from a group view}
# class GRE_LAYOUT Layout {use to access the layout being viewed if the script is run from a layout view}
# numeric ToolIsActive Will be 0 if tool is inactive or 1 if tool is active
#
# The following values are also predefined and are valid when the various On...()
# functions are called which deal with pointer and keyboard events.
# numeric PointerX Pointer X coordinate within view in pixels
# numeric PointerY Pointer Y coordinate within view in pixels
# numeric ShiftPressed 1 if <shift> key being pressed or 0 if not
# numeric CtrlPressed 1 if <ctrl> key being pressed or 0 if not
# numeric LeftButtonPressed 1 if left pointer button pressed or 0 if not
# numeric RightButtonPressed 1 if right pointer button pressed or 0 if not
# numeric MiddleButtonPressed 1 if middle pointer button pressed or 0 if not
#
##################################################################################
##################################################################################
#
# User Defined Global Variable
#
# Increase this variable value to expand the size of the Outcrop Trace Pattern that
# is created when using the tool.
# Example:
# Final Left Extents = Original Left - (width of box * stretchClipBoundry)
# Final Right Extents = Original Right + (width of box * stretchClipBoundry)
# Final Upper Extents = Original Upper + (height of box * stretchClipBoundry)
# Final Lower Extents = Original Lower - (height of box * stretchClipBoundry)
#
##################################################################################
numeric stretchClipBoundry = 0.5;
##################################################################################
##################################################################################
# ToolScript Global Variables
##################################################################################
class XMLDOC dlgdoc;
class XMLNODE node;
class GUI_DLG dlg;
class GUI_CTRL_LABEL promptString;
class GUI_CTRL_TOGGLEBUTTON overturnedToggle;
class GUI_CTRL_EDIT_NUMBER strikeFld, dipFld, dipDirFld;
class GUI_CTRL_PUSHBUTTON btnACCEPT, btnCANCEL, btnDELETE, btnCLOSE;
class GUI_FORM_RADIOGROUP rGpMODE;
class GUI_CTRL_LABEL aLINECOLOR_LABEL;
class GUI_CTRL_COLORBUTTON btnLINECOLOR;
class MdispPOLYLINETOOL tool;
class POLYLINE poly;
class RVC_RASTER DEM, PLANE, clipRast;
class VECTOR StrikeDipVector, OutcropPVector, intersectV, clipV;
class GEOREF SDgeoref, demGeoref;
class RVC_GEOREFERENCE rvcGeorefDEM;
class TRANSMODEL demCalibModel;
class TRANS2D_MAPGEN transDEMobjToMap;
class SR_COORDREFSYS demCRS, utmCRS;
class TRANS2D_MAPGEN transGeogUTM; # coordinate transformation from geographic to UTM coordinates
class GRE_GROUP activegroup, firstGroup;
class GRE_LAYER activeLayer;
class GRE_LAYER_RASTER rasterLayer, intersectRlayer;
class GRE_LAYER_VECTOR SDvectorLayer, OPvectorLayer, intersectVlayer, clipVlayer;
class GRE_VECTOR_POINTS SDpoints, OPpoints;
class GRE_VECTOR_LINES OPlines;
class GRE_VECTOR_POLYS OPpolys;
class DATABASE SDdb, OPlinedb, OPpointdb, iVpointdb, iVlinedb;
class DBTABLEINFO SDtable, OPlinetable, OPpointtable, iVpointtable, iVlinetable;
class DBFIELDINFO SDstationfield;
numeric demZscale; # z-scale value for DEM
numeric a=0, b=0, c=0, d=0; # linear coefficients for the equation of the plane
numeric dipX, dipY; # partial derivatives of dip (negative slope) in x and y directions
numeric dipAng; # magnitude of dip angle
numeric dipDir; # azimuth of dip direction (0 to 360 degrees clockwise from north)
numeric strikeAng; # azimuth of strike line
numeric station, overturned;
numeric setDefaultWhenClose;
numeric selectedElementNum;
numeric currentlyEditing = 0;
numeric waitingForAddConfirmation = 0;
numeric deletedStation;
numeric transGeogCoords; # flag to indicate map coordinates need to be transformed from
# geographic to UTM
array numeric records[1];
class POINT3D ctrPt;
string Current_Mode = "ADD";
string genericQuery = "(Outcrop.Station == ";
string currentQuery;
##################################################################################
##################################################################################
# Function Prototypes
##################################################################################
# Function prototypes with a class as a return value are unallowed at this time
##################################################################################
#func class POINT3D setXY(numeric num);
#func class POINT3D setZ(class POINT3D pt3d, numeric num);
#func class POINT3D computeDifference(class POINT3D pt1, class POINT3D pt2);
#func class POINT3D computeCenterPoint(class POINT3D pt1, class POINT3D pt2, class POINT3D pt3);
#func class POINT3D computePlanarCoefficients();
# Functions used to compute the Strike and Dip values
##################################################################################
proc createPlaneRaster(class POINT3D testPoint, class POINT3D point2, class POINT3D point3);
func Determinant(a1, b1, a2, b2);
func isZValid(numeric z);
proc computeVerticalPlane();
proc computeNonVerticalPlane();
proc computeStrikeDip();
# Functions used to manage the Strike and Dip data
##################################################################################
func addPointToVector(class POINT3D pt);
proc addRecordToDatabase( numeric pointElemNumAdded );
proc UpdateQuery( );
proc deleteStation();
# Functions used to initiliaze or maintain the display layers and databases
##################################################################################
func checkLayer();
proc checkGeoref();
proc createDestVectors();
func createVectorGeoreference();
func vectorLayerExists(string name);
proc addVectorsToDisplay();
proc initializeBeddingDatabaseTable();
proc createBeddingFields();
proc initializeOutcropDatabaseTable();
# Functions used to manage callbacks
##################################################################################
proc cbToolApply(class MdispPOLYLINETOOL tool);
proc cbLayer();
proc cbGroup();
proc cbClose();
# Functions called from interaction with the dialog
##################################################################################
proc OnModeChange( );
proc OnAccept( );
proc OnCancel( );
proc OnDelete();
proc OnClose();
proc OnOverturnedChanged();
# Basic ToolScript functions
##################################################################################
proc OnRightButtonPress();
proc OnLeftButtonPress();
proc OnInitialize ();
proc OnDestroy ();
proc OnActivate ();
proc OnDeactivate ();
##################################################################################
##################################################################################
# Function Definitions
##################################################################################
##################################################################################
#### Get the x and y screen coordinates from a polyline vertex
func class POINT3D setXY(numeric num) {
local class POINT3D tmp = poly.GetVertex(num - 1);
return tmp;
}
##################################################################################
#### Convert polyline vertex coordinates to DEM map coordinates. Get the z value
#### from the elevation raster at that location.
func class POINT3D setZ(class POINT3D pt3d) {
local class POINT2D tmp;
tmp.x = pt3d.x;
tmp.y = pt3d.y;
# convert from screen to view coords
tmp = TransPoint2D(tmp, ViewGetTransViewToScreen(View, 1));
# convert from view to map coords
tmp = TransPoint2D(tmp, ViewGetTransMapToView(View, rasterLayer.MapRegion.CoordRefSys, 1));
# convert from map coords to object
local class Georef georef;
local class POINT2D objcoord;
georef = GetLastUsedGeorefObject(DEM);
objcoord = MapToObject(georef, tmp.x, tmp.y, DEM);
# with the object coordinates we can get the cell value
numeric cellval = DEM[objcoord.y, objcoord.x];
# return a 3D point with map coordinates and elevation
local class POINT3D retPnt;
retPnt.x = tmp.x;
retPnt.y = tmp.y;
retPnt.z = cellval * demZscale; # a raster cell value with Z-scaling applied
return retPnt;
}
##################################################################################
#### Calculate the difference between two points
func class POINT3D computeDifference(class POINT3D pt1, class POINT3D pt2) {
local class POINT3D diff;
diff.x = pt1.x - pt2.x;
diff.y = pt1.y - pt2.y;
diff.z = pt1.z - pt2.z;
return diff;
}
##################################################################################
#### Compute the center point of the three points
func class POINT3D computeCenterPoint(class POINT3D pt1, class POINT3D pt2, class POINT3D pt3) {
local class POINT3D ctr;
ctr.x = (pt1.x + pt2.x + pt3.x) / 3;
ctr.y = (pt1.y + pt2.y + pt3.y) / 3;
ctr.z = (pt1.z + pt2.z + pt3.z) / 3;
return ctr;
}
##################################################################################
#### Convert geographic coordinates to planar coordinates in the designated UTM CRS
func class POINT3D convertGeogToUTM(class POINT3D geogPt) {
local class POINT2D mapCoords; # 2D map coordinates
local class POINT3D utmPt; # 3D point to return
mapCoords.x = geogPt.x; # read geographic coordinates from 3D point to 2D point
mapCoords.y = geogPt.y;
mapCoords = TransPoint2D(mapCoords, transGeogUTM); # convert 2D point to UTM
utmPt.x = mapCoords.x; # assign coordinate values to 3D point to return
utmPt.y = mapCoords.y;
utmPt.z = geogPt.z;
return utmPt;
}
##################################################################################
#### Convert planar coordinates in the designated UTM CRS back to geographic
func class POINT3D convertUTMtoGeog(class POINT3D utmPt) {
local class POINT2D mapCoords; # 2D map coordinates
local class POINT3D geogPt; # 3D point to return
mapCoords.x = utmPt.x;
mapCoords.y = utmPt.y;
mapCoords = TransPoint2D(mapCoords, transGeogUTM, 1); # inverse transformation
geogPt.x = mapCoords.x;
geogPt.y = mapCoords.y;
geogPt.z = utmPt.z;
return geogPt;
}
##################################################################################
#### Find angle between grid north and true north for station location (center point)
func findAngleToNorth (class POINT3D ctr) {
local class POINT2D mapCoords;
local class SR_COORDREFSYS ptCRS;
if (transGeogCoords == 1) {
ctr = convertGeogToUTM(ctrPt);
ptCRS = utmCRS;
}
else {
ptCRS = demCRS;
}
return ptCRS.ComputeAngleToNorth(ctr);
}
##################################################################################
#### Compute linear coefficients of the equation of the plane from the three
#### 3D points in the polyline tool
func class POINT3D computePlanarCoefficients( ) {
local class POINT3D pt1, pt2, pt3; # three points on a plane
local class POINT3D pt21, pt31; # differences between point positions
local class POINT3D ctrPt;
pt1 = setXY(1); # get screen coordinates for polyline vertex
pt1 = setZ(pt1); # convert x,y to map coordinates and set Z value from DEM
pt2 = setXY(2);
pt2 = setZ(pt2);
pt3 = setXY(3);
pt3 = setZ(pt3);
# if one of the z values is bad, then stop - return invalid ctrPt
if (!(isZValid(pt1.z) && isZValid(pt2.z) && isZValid(pt3.z)))
{
string s$ = "Error: One (or more) of the points has an invalid elevation.";
PopupMessage(s$);
ctrPt.x = NullValue(DEM);
ctrPt.y = NullValue(DEM);
ctrPt.z = NullValue(DEM);
return ctrPt;
}
else
{
# compute the center point of the triangle, this is the point that is created
# and attributed with the information about strike and dip.
ctrPt = computeCenterPoint(pt1, pt2, pt3);
# if DEM map coordinates are geographic, call function to convert each vertex point
# to planar UTM coordinates to compute equation of the plane
if (transGeogCoords == 1) {
pt1 = convertGeogToUTM(pt1);
pt2 = convertGeogToUTM(pt2);
pt3 = convertGeogToUTM(pt3);
}
# Compute linear coefficients for the equation of the plane
pt21 = computeDifference(pt2, pt1);
pt31 = computeDifference(pt3, pt1);
a = Determinant( pt21.y, pt21.z, pt31.y, pt31.z );
b = Determinant( pt21.x, pt21.z, pt31.x, pt31.z ) * -1;
c = Determinant( pt21.x, pt21.y, pt31.x, pt31.y );
# call function to compute outcrop trace
createPlaneRaster(pt1, pt2, pt3);
return ctrPt;
}
}
##################################################################################
#### Use the computed planar coefficients to determine the outcrop trace pattern
#### based on the plane intersecting the DEM raster
proc createPlaneRaster(class POINT3D testPoint, class POINT3D point2, class POINT3D point3) {
local numeric i, j, k, l, tempVal;
local class POINT2D tempPoint;
local numeric numLins, numCols;
local class RVC_RASTER RastPlane;
local class RVC_GEOREFERENCE georefPlane, georefClip;
local class POINT3D ulPt, urPt, lrPt, llPt;
local class CTRLPOINT ctrlPointList[4]; # array of control points
local class RVC_VECTOR tempVect, ClipVect;
# compute constant d (z-intercept) of plane using test point
d = (-1) * ((a*testPoint.x) + (b* testPoint.y) + (c*testPoint.z));
# Determine the extents of the polytool
local class RECT polyRect = poly.ComputeExtents();
local class TRANS2D_MAPGEN transScrnToLayer = View.GetTransLayerToScreen(rasterLayer, true);
ulPt = polyRect.pt1; # screen coordinates for upper left and lower right corners of polytool extents
lrPt = polyRect.pt2;
ulPt = transScrnToLayer.ConvertPoint2DFwd(ulPt); # extents in DEM object coordinates
lrPt = transScrnToLayer.ConvertPoint2DFwd(lrPt);
# Use the user-defined variable to expand the extents of the outcrop trace
numLins = (lrPt.y - ulPt.y) * stretchClipBoundry;
numCols = (lrPt.x - ulPt.x) * stretchClipBoundry;
ulPt.x = floor(ulPt.x - numCols);
ulPt.y = floor(ulPt.y - numLins);
lrPt.x = floor(lrPt.x + numCols);
lrPt.y = floor(lrPt.y + numLins);
# Make sure extents are within the bounds of the DEM raster
if (ulPt.x < 1) then ulPt.x = 1;
if (ulPt.y < 1) then ulPt.y = 1;
if (lrPt.x > DEM.$Info.NumCols) then lrPt.x = DEM.$Info.NumCols;
if (lrPt.y > DEM.$Info.NumLins) then lrPt.y = DEM.$Info.NumLins;
# Get number of DEM lines and columns in the outcrop trace area
numLins = lrPt.y - ulPt.y + 1;
numCols = lrPt.x - ulPt.x + 1;
# set raster object coordinates for upper right and lower left corners of trace area
urPt.x = lrPt.x; urPt.y = ulPt.y;
llPt.x = ulPt.x; llPt.y = lrPt.y;
# Create temporary binary raster for outcrop trace area with georefence using same CRS as DEM
CreateTempRaster(RastPlane, numLins, numCols, "binary");
georefPlane.SetCoordRefSys(demCRS);
georefPlane.SetCalibModel(demCalibModel);
# create georeference control points for plane raster
ctrlPointList[1].source.x = 1; # upper left control point
ctrlPointList[1].source.y = 1;
ctrlPointList[1].dest = transDEMobjToMap.ConvertPoint2DFwd(ulPt);
ctrlPointList[2].source.x = numCols; # upper right control point
ctrlPointList[2].source.y = 1;
ctrlPointList[2].dest = transDEMobjToMap.ConvertPoint2DFwd(urPt);
ctrlPointList[3].source.x = numCols; # lower right control point
ctrlPointList[3].source.y = numLins;
ctrlPointList[3].dest = transDEMobjToMap.ConvertPoint2DFwd(lrPt);
ctrlPointList[4].source.x = 1; # lower left control point
ctrlPointList[4].source.y = numLins;
ctrlPointList[4].dest = transDEMobjToMap.ConvertPoint2DFwd(llPt);
# set control point list for georeference and assign georeference to RastPlane
georefPlane.SetPointList(ctrlPointList);
georefPlane.Make(RastPlane);
# Set up a status message to inform user of progress
local class STATUSDIALOG statusD;
statusD.Create();
local class STATUSCONTEXT statusC = statusD.CreateContext();
StatusSetMessage(statusC, "Computing Strike and Dip Values and Outcrop Pattern");
StatusSetBar(statusC, i, DEM.$Info.NumLins);
# Loop through cells in RastPlane; get Z value for corresponding cell in DEM and
# use its map coordinates to compute Z value from the equation for the strike/dip plane.
# If Plane is above ground, set cell value in RastPlane to 1. After all cells are
# processed, boundary between 1's and 0's is the outcrop trace for the plane.
for i = 1 to RastPlane.$Info.NumLins {
StatusSetBar(statusC, i, DEM.$Info.NumLins);
for j = 1 to RastPlane.$Info.NumCols {
k = ulPt.y + i - 1; l = ulPt.x + j - 1; # matching object coordinates in DEM
tempPoint.x = l;
tempPoint.y = k;
tempPoint = transDEMobjToMap.ConvertPoint2DFwd(tempPoint); # cell location in DEM map coordinates
if (transGeogCoords == 1) { # if DEM has geographic CRS, convert map coords to UTM to match plane coordinates
tempPoint = TransPoint2D(tempPoint, transGeogUTM);
}
tempVal = (-1) * (a * tempPoint.x + b * tempPoint.y + d) / c;
if (tempVal > DEM[k,l] * demZscale) then
RastPlane[i, j] = 1;
}
}
# Create a Vector containing the outcrop trace line. This is created by finding the boundaries
# between the 0's and 1's in the raster that compared the DEM to the computed plane.
CreateTempVector(tempVect);
tempVect = RasterToVectorBound(RastPlane, "DoImpliedGeoref");
# Create clip area to extract outcrop trace line and leave border behind;
# set corners of clip area in map coordinates, equivalent to 1 DEM cell smaller than trace area
local class POINT2D rastScale;
DEM.ComputeScaleFromGeoref(rvcGeorefDEM, rastScale, 0); # x and y map dimensions of a DEM cell
local class POINT2D vertexPt;
local class POLYLINE polyL;
vertexPt.x = ctrlPointList[1].dest.x; # upper left vertex
vertexPt.y = ctrlPointList[1].dest.y;
polyL.AppendVertex(vertexPt);
vertexPt.x = ctrlPointList[2].dest.x - rastScale.x; # upper right vertex
vertexPt.y = ctrlPointList[2].dest.y;
polyL.AppendVertex(vertexPt);
vertexPt.x = ctrlPointList[3].dest.x - rastScale.x; # lower right vertex
vertexPt.y = ctrlPointList[3].dest.y + rastScale.y;
polyL.AppendVertex(vertexPt);
vertexPt.x = ctrlPointList[4].dest.x; # lower left vertex
vertexPt.y = ctrlPointList[4].dest.y + rastScale.y;
polyL.AppendVertex(vertexPt);
polyL.SetClosed(1); # close the polyline
polyRect = polyL.ComputeExtents();
CreateTempVector(ClipVect, "VectorTookit,Polygonal","NoDatabase", polyRect);
local class RVC_GEOREFERENCE georefClipVect;
georefClipVect.SetCoordRefSys(demCRS);
georefClipVect.SetImplied();
georefClipVect.Make(ClipVect);
VectorAddPolyLine(ClipVect, polyL);
VectorValidate(ClipVect);
# Extract outcrop trace line to new vector using the ClipVect to leave behind the border created
# by the autobounds procedure.
CreateTempVector(intersectV);
intersectV = VectorExtract(ClipVect, tempVect, "InsideClip", "", "1", "1", "0");
# if DEM has geographic coordinates, reproject UTM points passed to this function back to geographic
# so they can be added to the outcrop trace vector
if (transGeogCoords == 1) {
testPoint = convertUTMtoGeog(testPoint);
point2 = convertUTMtoGeog(point2);
point3 = convertUTMtoGeog(point3);
}
# Add the three points selected by the user
VectorAddPoint(intersectV, testPoint.x, testPoint.y, testPoint.z);
VectorAddPoint(intersectV, point2.x, point2.y, point2.z);
VectorAddPoint(intersectV, point3.x, point3.y, point3.z);
# Create a Point Database table and record with the station number and color values
iVpointdb = OpenVectorPointDatabase(intersectV);
iVpointtable = TableCreate(iVpointdb, "Outcrop", "Intersecting Outcrop Pattern");
TableAddFieldInteger(iVpointtable, "Station", 4);
TableAddFieldInteger(iVpointtable, "RedComp", 3);
TableAddFieldInteger(iVpointtable, "GreenComp", 3);
TableAddFieldInteger(iVpointtable, "BlueComp", 3);
local numeric recNumPoints = TableWriteRecord(iVpointtable, 0, -999, 0, 0, 0);
# Create a Line Database table and with the station number and color values
iVlinedb = OpenVectorLineDatabase(intersectV);
iVlinetable = TableCreate(iVlinedb, "Outcrop", "Intersecting Outcrop Pattern");
TableAddFieldInteger(iVlinetable, "Station", 4);
TableAddFieldInteger(iVlinetable, "RedComp", 3);
TableAddFieldInteger(iVlinetable, "GreenComp", 3);
TableAddFieldInteger(iVlinetable, "BlueComp", 3);
local numeric recNumLines = TableWriteRecord(iVlinetable, 0, -999, 0, 0, 0);
# Attach all the points to the station
local numeric numPoints = NumVectorPoints(intersectV);
local numeric k;
local array numeric rec[1];
rec[1] = recNumPoints;
for k = 1 to numPoints {
TableWriteAttachment(iVpointtable, k, rec, 1, "point");
}
# Attach all the lines to the station
local numeric numLines = NumVectorLines(intersectV);
rec[1] = recNumLines;
for k = 1 to numLines {
TableWriteAttachment(iVlinetable, k, rec, 1, "line");
}
# Add the outcrop trace line to the view for the user to evaluate
intersectVlayer = GroupQuickAddVectorVar(firstGroup, intersectV);
intersectVlayer.Point.NormalStyle.DrawCircle = 1;
intersectVlayer.Point.NormalStyle.Fill = 1;
intersectVlayer.Line.StyleMode = "AllSame";
intersectVlayer.Line.NormalStyle.Color = btnLINECOLOR.GetColor();
intersectVlayer.Line.NormalStyle.Width = 1.05;
ViewRedrawLayer(View, intersectVlayer);
# Clean up the resources used
statusC.Clear();
statusD.Destroy();
CloseVector(tempVect);
CloseVector(ClipVect);
CloseRaster(RastPlane);
# CloseRaster(clipRast);
# CloseRaster(PLANE);
}
##################################################################################
#### Function to return the determinant of a 2 x 2 matrix formed by the four
#### numbers in the argument. The numbers can be in either row or column order;
#### the result is the same in either order.
func Determinant(a1, b1, a2, b2) {
local det;
det = a1 * b2 - a2 * b1 ;
return det;
}
##################################################################################
#### Check that the z value is not null (either null or beyond rast extents)
func isZValid(numeric z) {
local numeric nullval = NullValue(DEM);
return z != nullval;
}
##################################################################################
#### Compute the strike and dip for a vertical plane
proc computeVerticalPlane( ) {
local numeric northAng;
northAng = findAngleToNorth(ctrPt);
dipAng = 90;
if ( a == 0 )
{
# If a = b = c = 0, points are colinear, can't find plane;
# tool script should trap this as an error and alert the user
# to repick points -- tool shouldn't let this happen -djg
if ( b == 0 )
{
strikeAng = null;
dipDir = null;
dipAng = null;
PopupMessage("Points are colinear and so don't define a unique plane.");
}
# Plane is also parallel to x-axis, raw strike = 90 and arctan(-b/a) is undefined
else
{
strikeAng = 90 + northAng;
}
}
else
{
# Raw strike angle = arctan of the partial derivative of x versus y = (-b / a)
# (this is slope of line of intersection of plane with horizontal = strike line
# in plot of x vs y, which gives the angle equivalent to azimuth in conventional
# y vs x plot)
strikeAng = atand( -b / a ) + northAng;
if ( strikeAng < 0 ) then
strikeAng = 180 + strikeAng;
}
}
##################################################################################
#### Compute the strike and dip for a nonvertical plane
proc computeNonVerticalPlane( ) {
local numeric northAng;
northAng = findAngleToNorth(ctrPt);
# Partial derivative of dip in x-direction (opposite the slope component)
dipX = a / c;
# Plane is parallel to y-axis (north-south), dipY is 0;
if ( b == 0 ) {
# Plane is parallel to both x and y and thus is horizontal
if ( a == 0 ) {
dipAng = 0;
strikeAng = null; # strikeAng is undefined
dipDir = null;
# Plane is not horizontal, raw strike is north-south
} else {
dipAng = atand( dipX );
if ( dipAng < 0 ) {
dipAng = -1 * dipAng;
}
# Dip direction is due east
if ( dipX > 0 ) {
dipDir = 90 + northAng;
strikeAng = 0 + northAng;
} else {
dipDir = 270 + northAng;
strikeAng = 180 + northAng;
}
}
# Strike not parallel to y-axis, dipY is nonzero
} else {
# Partial derivative of dip in y-direction (opposite the slope component)
dipY = b / c;
# Compute dip angle from dipX and dipY
dipAng = atand( sqrt( sqr(dipX) + sqr(dipY) ) );
if ( dipAng < 0 ) {
dipAng = -1 * dipAng;
}
# Compute azimuth of dip direction from dipX and dipY; value returned by atand() is between -90 and +90,
# so must adjust result to get positive azimuth value between 0 and 360 depending on which quadrant dip
# vector is in.
dipDir = atand( dipX / dipY ) + northAng;
# Dip direction is due north or south
if ( dipDir == 0 ) {
# Dip direction is south
if ( dipY < 0 ) {
dipDir = 180;
}
}
# Dip direction in SE or NW quadrants, must adjust value
if ( dipDir < 0 ) {
# Dip direction in SE quadrant
if ( dipX > 0 ) {
dipDir = 180 + dipDir;
# Dip direction in NW quadrant
} else {
dipDir = 360 + dipDir;
}
# Dip direction in SW quadrant, must adjust value
} else if ( dipX < 0 ) {
dipDir = 180 + dipDir;
}
# Compute strike angle from dip direction and convert negative values to
# positive azimuth in range 0 to 360
# Strike direction from right-hand rule
strikeAng = dipDir - 90;
if ( strikeAng < 0 ) {
strikeAng = 360 + strikeAng;
}
}
}
##################################################################################
#### Compute strike and dip values based on the user input
proc computeStrikeDip( ) {
local numeric isVertical = (c==0);
if (isVertical) {
computeVerticalPlane();
} else {
computeNonVerticalPlane();
}
}
##################################################################################
#### Add the point to the StrikeDipVector and return its element number
func addPointToVector(class POINT3D pt) {
local numeric x = pt.x;
local numeric y = pt.y;
local numeric z = pt.z;
VectorAddPoint(StrikeDipVector, x, y, z);
VectorValidate(StrikeDipVector);
return FindClosestPoint(StrikeDipVector, x, y, SDgeoref);
}
##################################################################################
#### Add a record to the StrikeDipVector point database with appropriate values
proc addRecordToDatabase( numeric pointElemNumAdded ) {
local class POINT3D pt1, pt2, pt3;
# Get the original three points from the polytool
pt1 = setXY(1);
pt1 = setZ(pt1);
pt2 = setXY(2);
pt2 = setZ(pt2);
pt3 = setXY(3);
pt3 = setZ(pt3);
# If editing, use the same station number
if (Current_Mode == "EDIT") {
station = deletedStation;
# If adding a new one, set as the highest station incremented by 1
} else {
station = 0;
local numeric tempVal;
local numeric numRecords = NumRecords(SDtable);
local numeric i;
for i = 1 to numRecords {
tempVal = TableReadFieldNum(SDtable, "Station", i);
if (tempVal > station) {
station = tempVal;
}
}
station = station + 1;
}
# Round and save Strike Dip Values as Integers
strikeAng = round(strikeAng);
dipAng = round(dipAng);
dipDir = round(dipDir);
# Write the record for the Strike Dip point and attach it
records[1] = TableWriteRecord(SDtable, 0, station, strikeAng, dipAng, dipDir, overturned,
pt1.x, pt1.y, pt2.x, pt2.y, pt3.x, pt3.y, 0);
TableWriteAttachment(SDtable, pointElemNumAdded, records, 1, "point");
# Write records into the temporary outcrop trace pattern vector. They will be transfered
# to the actual outcrop trace vector in the next step.
local class COLOR tempColor = btnLINECOLOR.GetColor();
TableWriteRecord(iVpointtable, 1, station, tempColor.red, tempColor.green, tempColor.blue);
TableWriteRecord(iVlinetable, 1, station, tempColor.red, tempColor.green, tempColor.blue);
# Vector Copy Elements doesn't like copying to an object that already has elements, so if
# it does, we introduce a temporary vector before using the VectMerge process.
if ( NumVectorPoints(OutcropPVector) + NumVectorLines(OutcropPVector) > 0 ) {
# Close the databases and remove the layers
CloseDatabase(SDdb);
CloseDatabase(OPpointdb);
CloseDatabase(OPlinedb);
GroupRemoveLayer(firstGroup, SDvectorLayer);
GroupRemoveLayer(firstGroup, OPvectorLayer);
# Copy the temporary outcrop trace vector to a temporary vector, then merge it with the
# final destination outcrop trace pattern vector
local VECTOR tempMergeVector;
CreateTempVector(tempMergeVector);
VectorCopyElements(OutcropPVector, tempMergeVector);
OutcropPVector = VectMerge(tempMergeVector, intersectV);
# Add the layers back to the display;
addVectorsToDisplay();
# Reinitialize and clean up the vector databases
initializeBeddingDatabaseTable();
initializeOutcropDatabaseTable();
TableRemoveDuplicateRecords(OPpointtable);
TableRemoveDuplicateRecords(OPlinetable);
CloseVector(tempMergeVector);
# If the outcrop trace vector has no elements, just copy straight to it.
} else {
# Close the databases and remove the layers
CloseDatabase(SDdb);
CloseDatabase(OPpointdb);
CloseDatabase(OPlinedb);
GroupRemoveLayer(firstGroup, SDvectorLayer);
GroupRemoveLayer(firstGroup, OPvectorLayer);
VectorCopyElements(intersectV, OutcropPVector);
# Add the layers back to the display;
addVectorsToDisplay();
# Reinitialize and clean up the vector databases
initializeBeddingDatabaseTable();
initializeOutcropDatabaseTable();
TableRemoveDuplicateRecords(OPpointtable);
TableRemoveDuplicateRecords(OPlinetable);
}
}
##################################################################################
#### This function creates the Query used to display the appropriate lines attached
#### to specific stations. It also manages the color they appear in.
proc UpdateQuery( ) {
currentQuery = genericQuery + "-1)";
# Scan through every record in the StrikeDipVector Point database
local numeric numPoints = NumVectorPoints(StrikeDipVector);
local array numeric records[10];
local numeric i;
for i = 1 to numPoints {
TableReadAttachment(SDtable, i, records, "point");
# If the record value states that it should be displayed, add an extra OR statement
# to the current query so it is displayed.
if (TableReadFieldNum(SDtable, "DisplayOPdata", records[1]) != 0) {
currentQuery = currentQuery + " || " + genericQuery
+ NumToStr(TableReadFieldNum(SDtable, "Station", records[1])) + ")";
}
}
# Set the current query and redraw the layer
OPpoints.Select.Query = currentQuery;
OPlines.Select.Query = currentQuery;
ViewRedrawLayer(View, OPvectorLayer);
}
##################################################################################
#### This function is designed to delete all lines, points and database records that
#### are related to the station that is being deleted.
proc deleteStation( ) {
local numeric numToDelete = 0;
local numeric i;
local array numeric tempRecords[10];
local numeric tempStation;
local numeric numRecords;
# Turn off the highlighted point
SDpoints.HighlightSingle(selectedElementNum, "Toggle");
# Get the station number from the selected point
deletedStation = StrikeDipVector.point[selectedElementNum].Bedding.Station;
# Scan through and mark all lines attached to the deleted station to be deleted.
local numeric numLines = NumVectorLines(OutcropPVector);
local array numeric linesToDelete[numLines];
for i = 1 to numLines {
tempStation = OutcropPVector.line[i].Outcrop.Station;
if (tempStation == deletedStation) {
# Unattach the records if the line is being marked to delete
numRecords = TableReadAttachment(OPlinetable, i, tempRecords, "line");
TableRemoveAttachment(OPlinetable, i, tempRecords, numRecords, "line");
# Save line number to delete
numToDelete = numToDelete + 1;
linesToDelete[numToDelete] = i;
}
}
# Delete all the marked lines
if (numToDelete > 0) {
VectorDeleteLines(OutcropPVector, linesToDelete, numToDelete);
}
numToDelete = 0;
# Scan through and mark all points attached to the deleted station to be deleted.
local numeric numPoints = NumVectorPoints(OutcropPVector);
local array numeric pointsToDelete[numPoints];
for i = 1 to numPoints {
tempStation = OutcropPVector.point[i].Outcrop.Station;
if (tempStation == deletedStation) {
# Unattach the records if the point is being marked to delete
numRecords = TableReadAttachment(OPpointtable, i, tempRecords, "point");
TableRemoveAttachment(OPpointtable, i, tempRecords, numRecords, "point");
# Save point number to delete
numToDelete = numToDelete + 1;
pointsToDelete[numToDelete] = i;
}
}
# Delete all the marked points
if (numToDelete > 0) {
VectorDeletePoints(OutcropPVector, pointsToDelete, numToDelete);
}
numToDelete = 0;
# Scan through and delete the Strike and Dip station
numPoints = NumVectorPoints(StrikeDipVector);
local numeric pointToDelete;
for i = 1 to numPoints {
tempStation = StrikeDipVector.point[i].Bedding.Station;
if (tempStation == deletedStation) {
# Unnattach the records if the point is being marked to delete
numRecords = TableReadAttachment(SDtable, i, tempRecords, "point");
TableRemoveAttachment(SDtable, i, tempRecords, numRecords, "point");
# Save point number to delete
numToDelete = numToDelete + 1;
pointsToDelete[numToDelete] = i;
}
}
# Delete all the marked points
if (numToDelete > 0) {
VectorDeletePoints(StrikeDipVector, pointsToDelete, numToDelete);
}
# Remove all the records that are no longer attached
TableRemoveUnattachedRecords(OPlinetable);
TableRemoveUnattachedRecords(OPpointtable);
TableRemoveUnattachedRecords(SDtable);
VectorValidate(OutcropPVector);
VectorValidate(StrikeDipVector);
}
##################################################################################
#### Checks the first layer to see if it is valid DEM raster
func checkLayer( ) {
local numeric valid = false;
local class GRE_LAYER layer = firstGroup.FirstLayer;
if (layer.TypeID == "Surface") layer = layer.NextLayer;
# Get name of active layer if it is usable. If not output an error message.
if (layer.TypeID != "")
{
if (layer.TypeID == "Raster")
{
DispGetRasterFromLayer(DEM, layer);
if (DEM.$Info.Type == "binary" or
DEM.$Info.Type == "4-bit unsigned" or
DEM.$Info.Type == "8-bit signed" or
DEM.$Info.Type == "8-bit unsigned" or
DEM.$Info.Type == "16-bit signed" or
DEM.$Info.Type == "16-bit unsigned" or
DEM.$Info.Type == "32-bit signed" or
DEM.$Info.Type == "32-bit unsigned" or
DEM.$Info.Type == "32-bit floating-point" or
DEM.$Info.Type == "64-bit signed" or
DEM.$Info.Type == "64-bit unsigned" or
DEM.$Info.Type == "64-bit floating-point")
{
rasterLayer = layer;
valid = true;
}
}
}
return valid;
}
##################################################################################
#### Creates the vector objects to store the resultant points and outcrop trace lines
proc createDestVectors( ) {
local string flags$ = "Polygonal,3DVector";
GetOutputVector(StrikeDipVector, flags$);
GetOutputVector(OutcropPVector, flags$);
createVectorGeoreference();
}
##################################################################################
#### Creates the vector georeference from the raster if necessary
func createVectorGeoreference( ) {
local class Georef rgeoref;
# determine if any georef objects exist
local string file$ = StrikeDipVector.$INFO.Filename;
local string type$ = "GEOREF";
local numeric parent = GetObjectNumber(StrikeDipVector);
local class STRINGLIST sl = GetAllObjectNames(file$, type$, parent);
local numeric hasGeoref = (sl.GetNumItems()!=0);
# If it has a Georef, use it. If it doesn't, create one from the DEM
if (hasGeoref) {
SDgeoref = GetLastUsedGeorefObject(StrikeDipVector);
} else {
rgeoref = GetLastUsedGeorefObject(DEM);
CreateImpliedGeoref(StrikeDipVector, rgeoref.CoordRefSys);
SDgeoref = GetLastUsedGeorefObject(StrikeDipVector);
CreateImpliedGeoref(OutcropPVector, rgeoref.CoordRefSys);
}
}
##################################################################################
#### This function is used to check if the vector is already in a layer in the View
func vectorLayerExists( string name ) {
local class VECTOR V;
# Set the active layer to the top layer in the first group
firstGroup.SetActiveLayer(firstGroup.LastLayer);
# Cycle through all the layers in the group
while (firstGroup.ActiveLayer != 0) {
if (firstGroup.ActiveLayer.TypeID == "Vector") {
DispGetVectorFromLayer(V, firstGroup.ActiveLayer);
# If the vector layer has the same name, return true
if (V.$Info.Name == name) {
CloseVector(V);
return true;
}
CloseVector(V);
}
# If there are no more layers to check, return false
if (firstGroup.ActiveLayer.PrevLayer == 0) {
return false;
}
# Set the active layer to the next one down
firstGroup.SetActiveLayer(firstGroup.ActiveLayer.PrevLayer);
}
return false;
}
##################################################################################
#### Add the vector to the display group if it is not already there
proc addVectorsToDisplay( ) {
# Check if the StrikeDipVector layer does not already exist
if(!vectorLayerExists(StrikeDipVector.$Info.Name))
{
# Add the vector to the View
SDvectorLayer = GroupQuickAddVectorVar(firstGroup, StrikeDipVector);
SDvectorLayer.IgnoreExtents = 1;
} else {
SDvectorLayer = firstGroup.ActiveLayer;
CloseVector(StrikeDipVector);
DispGetVectorFromLayer(StrikeDipVector, SDvectorLayer);
SDvectorLayer.IgnoreExtents = 1;
}
# Check if the OutcropPVector layer does not already exist
if(!vectorLayerExists(OutcropPVector.$Info.Name))
{
# Add the vector to the View
OPvectorLayer = GroupQuickAddVectorVar(firstGroup, OutcropPVector);
OPvectorLayer.IgnoreExtents = 1;
} else {
OPvectorLayer = firstGroup.ActiveLayer;
CloseVector(OutcropPVector);
DispGetVectorFromLayer(OutcropPVector, OPvectorLayer);
OPvectorLayer.IgnoreExtents = 1;
}
# Set up how the layers display points, lines and polygons
SDpoints = SDvectorLayer.Point;
OPpoints = OPvectorLayer.Point;
OPlines = OPvectorLayer.Line;
OPpolys = OPvectorLayer.Poly;
# Use the modified bedding.qry CartoScript developed by MicroImages to
# be used by this ToolScript if it exists in the same directory
SDpoints.Select.Mode = "All";
if (fexists(_context.ScriptDir + "/beddingStrikeDip.qry")) {
SDvectorLayer.Point.StyleMode = "ByScript";
string script = ScriptResourceReadFull("beddingStrikeDip.qry");
SDvectorLayer.Point.Script = script;
# If it doesn't exist, use large dots
} else {
SDvectorLayer.Point.StyleMode = "AllSame";
SDvectorLayer.Point.NormalStyle.xsize = 4;
SDvectorLayer.Point.NormalStyle.ysize = 4;
SDvectorLayer.Point.NormalStyle.DrawCircle = 1;
SDvectorLayer.Point.NormalStyle.Fill = 1;
}
# Create a script that will style the points by color
local string script = "DrawColor$ = NumToStr(Outcrop.RedComp) + \" \" + NumToStr(Outcrop.GreenComp)"
+ " + \" \" + NumToStr(Outcrop.BlueComp);\n";
# Set up the points be selected by query and styled by script
OPpoints.Select.Mode = "ByQuery";
OPpoints.Select.Query = "Outcrop.Station == -1";
OPpoints.StyleMode = "ByScript";
OPpoints.Script = script;
# Add to the script to make the lines thicker
script = script + "LineScale = 1;";
# Set up the lines be selected by query and styled by script
OPlines.Select.Mode = "ByQuery";
OPlines.Select.Query = "Outcrop.Station == -1";
OPlines.StyleMode = "ByScript";
OPlines.Script = script;
# Turn off polygons
OPpolys.Select.Mode = "None";
}
##################################################################################
#### Create the vector point database tables for the StrikeDipVector
proc initializeBeddingDatabaseTable( ) {
local string name$ = "Bedding";
local string desc$ = "Strike and dip bedding";
SDdb = OpenVectorPointDatabase(StrikeDipVector);
# If the table doesn't exist it needs to be created
if( !TableExists(SDdb, name$) ) {
SDtable = TableCreate(SDdb, name$, desc$);
SDtable.OneRecordPerElement = 1;
# Create the fields
createBeddingFields();
# Else the table does exist, get the info for it
} else {
SDtable = DatabaseGetTableInfo(SDdb, name$);
}
SDstationfield = FieldGetInfoByName(SDtable, "Station");
}
##################################################################################
#### Create the fields in the bedding table
proc createBeddingFields( ) {
TableAddFieldInteger(SDtable, "Station", 4);
TableAddFieldInteger(SDtable, "StrikeAngle", 4);
TableAddFieldInteger(SDtable, "DipAngle", 4);
TableAddFieldInteger(SDtable, "DipDirection", 4);
TableAddFieldInteger(SDtable, "Overturned", 4);
TableAddFieldFloat(SDtable, "ToolPos1_x", 4, 4);
TableAddFieldFloat(SDtable, "ToolPos1_y", 4, 4);
TableAddFieldFloat(SDtable, "ToolPos2_x", 4, 4);
TableAddFieldFloat(SDtable, "ToolPos2_y", 4, 4);
TableAddFieldFloat(SDtable, "ToolPos3_x", 4, 4);
TableAddFieldFloat(SDtable, "ToolPos3_y", 4, 4);
TableAddFieldInteger(SDtable, "DisplayOPdata", 2);
}
##################################################################################
#### Create the vector point database tables for the OutcropPVector
proc initializeOutcropDatabaseTable( ) {
local string name$ = "Outcrop";
local string desc$ = "Intersecting Outcrop Pattern";
OPpointdb = OpenVectorPointDatabase(OutcropPVector);
OPlinedb = OpenVectorLineDatabase(OutcropPVector);
# If the table doesn't exist it needs to be created
if( !TableExists(OPpointdb, name$) ) {
OPpointtable = TableCreate(OPpointdb, name$, desc$);
OPpointtable.OneRecordPerElement = 1;
# Add the fields in the table
TableAddFieldInteger(OPpointtable, "Station", 4);
TableAddFieldInteger(OPpointtable, "RedComp", 3);
TableAddFieldInteger(OPpointtable, "GreenComp", 3);
TableAddFieldInteger(OPpointtable, "BlueComp", 3);
# Else the table does exist, get the info for it
} else {
OPpointtable = DatabaseGetTableInfo(OPpointdb, name$);
}
# If the table doesn't exist it needs to be created
if( !TableExists(OPlinedb, name$) ) {
OPlinetable = TableCreate(OPlinedb, name$, desc$);
OPlinetable.OneRecordPerElement = 1;
# Add the fields in the table
TableAddFieldInteger(OPlinetable, "Station", 4);
TableAddFieldInteger(OPlinetable, "RedComp", 3);
TableAddFieldInteger(OPlinetable, "GreenComp", 3);
TableAddFieldInteger(OPlinetable, "BlueComp", 3);
# else the table does exist, get info for it
} else {
OPlinetable = DatabaseGetTableInfo(OPlinedb, name$);
}
}
##################################################################################
#### Callback for Apply button press for polyline (right click as well).
proc cbToolApply(class MdispPOLYLINETOOL tool)
{
# Make sure the bottom layer is still the DEM
if ( checkLayer() ) {
promptString.SetLabel("Computing Data...");
poly = tool.GetPolygon();
numeric numPoints = poly.GetNumPoints();
# If there are not exactly 3 points in accepted polyline
if ( numPoints != 3 ) {
string message$ = "Polyline must have exactly 3 points: "+ NumToStr(numPoints) + " found.";
PopupMessage(message$);
# Else there are 3 good points so do computations
} else {
# If the three points have been modified during the add process and are being
# recomputed, remove the original results
if (waitingForAddConfirmation) {
GroupRemoveLayer(firstGroup, intersectVlayer);
CloseVector(intersectV);
}
# Compute the coefficients from the 3 points
ctrPt = computePlanarCoefficients();
if (isZValid(ctrPt.z)) {
computeStrikeDip();
waitingForAddConfirmation = 1;
# Display the resulting Strike and Dip values in the dialog
strikeAng = round(strikeAng);
strikeFld.SetValue(strikeAng, 0);
dipFld.SetValue(round(dipAng), 0);
dipDirFld.SetValue(round(dipDir), 0);
overturnedToggle.SetEnabled(1);
# Manage the controls in the display
if (Current_Mode == "ADD") {
overturnedToggle.SetValue(0,0);
}
aLINECOLOR_LABEL.SetEnabled(1);
btnLINECOLOR.SetEnabled(1);
btnACCEPT.SetEnabled(1);
btnCANCEL.SetEnabled(1);
if (Current_Mode == "ADD") {
rGpMODE.SetEnabled(0);
btnCLOSE.SetEnabled(0);
promptString.SetLabel("Accept this New point?");
} else if (Current_Mode == "EDIT") {
promptString.SetLabel("Accept this Edited point?");
rGpMODE.SetEnabled(0);
btnCLOSE.SetEnabled(0);
}
}
}
# Else the bottom layer is not a valid DEM raster
} else {
string message$ = "The First Layer: "+activeLayer.Name+" is invalid.";
PopupMessage(message$);
}
}
##################################################################################
#### Callback for when the active layer changes
proc cbLayer( ) {
checkLayer();
checkGeoref();
}
##################################################################################
#### Callback for when the active group changes
proc cbGroup( ) {
activegroup = Layout.ActiveGroup;
WidgetAddCallback(activegroup.LayerSelectedCallback, cbLayer);
cbLayer();
}
##################################################################################
#### Close the window, switching to default tool
proc cbClose( ) {
tool.Managed = 0;
if (setDefaultWhenClose) {
setDefaultWhenClose = false;
View.SetDefaultTool();
}
}
##################################################################################
#### Manages the control dialog when the active mode has been changed
proc OnModeChange( ) {
# Get the current Mode
Current_Mode = rGpMODE.GetSelected();
# Set controls to default values
strikeFld.SetValue(0,0);
dipFld.SetValue(0, 0);
dipDirFld.SetValue(0, 0);
overturnedToggle.SetEnabled(1);
overturnedToggle.SetValue(0,0);
overturnedToggle.SetEnabled(0);
aLINECOLOR_LABEL.SetEnabled(0);
btnLINECOLOR.SetEnabled(0);
btnACCEPT.SetEnabled(0);
btnCANCEL.SetEnabled(0);
btnDELETE.SetEnabled(0);
# Based on Mode, set the controls and tool to appropriate settings
if ( Current_Mode == "VIEW") {
promptString.SetLabel("Select point to Toggle Data ON/OFF");
tool.managed = 0;
} else if ( Current_Mode == "ADD" ) {
aLINECOLOR_LABEL.SetEnabled(1);
btnLINECOLOR.SetEnabled(1);
promptString.SetLabel("Select 3 new points");
poly.Clear();
tool.SetPolygon(poly);
tool.managed = 1;
} else if ( Current_Mode == "EDIT" ) {
promptString.SetLabel("Select Point to Edit");
tool.managed = 0;
} else if ( Current_Mode == "DELETE" ) {
promptString.SetLabel("Select Point to Delete");
tool.managed = 0;
}
}
##################################################################################
#### Called when the user clicks on the Accept button in the dialog. Actions are
#### based on the current Mode of the ToolScript
proc OnAccept( ) {
# If in VIEW mode, toggle the display values of the selected station and update
if (Current_Mode == "VIEW") {
promptString.SetLabel("Resetting Query...");
btnACCEPT.SetEnabled(0);
btnCANCEL.SetEnabled(0);
# Toggle the values in the database
local array numeric records[10];
TableReadAttachment(SDtable, selectedElementNum, records, "point");
if (TableReadFieldNum(SDtable, "DisplayOPdata", records[1]) == 0) {
TableWriteField(SDtable, records[1], "DisplayOPdata", 1);
} else {
TableWriteField(SDtable, records[1], "DisplayOPdata", 0);
}
UpdateQuery();
SDpoints.HighlightSingle(selectedElementNum, "Toggle");
selectedElementNum = 0;
rGpMODE.SetEnabled(1);
btnCLOSE.SetEnabled(1);
promptString.SetLabel("Select point to Toggle Data ON/OFF");
# If in ADD Mode, add the computed Strike and Dip data to the Destination Vectors
} else if (Current_Mode == "ADD") {
promptString.SetLabel("Saving Data...");
overturned = overturnedToggle.GetValue();
overturnedToggle.SetEnabled(0);
aLINECOLOR_LABEL.SetEnabled(0);
btnLINECOLOR.SetEnabled(0);
btnACCEPT.SetEnabled(0);
btnCANCEL.SetEnabled(0);
# Remove the temporary vectors and save the Strike and Dip data
GroupRemoveLayer(firstGroup, intersectVlayer);
addRecordToDatabase( addPointToVector(ctrPt) );
CloseVector(intersectV);
waitingForAddConfirmation = 0;
# Reset the polytool
poly.Clear();
tool.SetPolygon(poly);
tool.Managed = 1;
# Update the display
UpdateQuery();
ViewRedrawLayer(View, SDvectorLayer);
selectedElementNum = 0;
aLINECOLOR_LABEL.SetEnabled(1);
btnLINECOLOR.SetEnabled(1);
rGpMODE.SetEnabled(1);
btnCLOSE.SetEnabled(1);
strikeFld.SetValue(0, 0);
dipFld.SetValue(0, 0);
dipDirFld.SetValue(0, 0);
promptString.SetLabel("Select 3 new points");
# If in EDIT Mode, reload the data or save the new edited data
} else if (Current_Mode == "EDIT") {
# Check if the point has already been selected and in modifying mode
if ( currentlyEditing ) {
promptString.SetLabel("Editing Data...");
overturned = overturnedToggle.GetValue();
overturnedToggle.SetEnabled(0);
aLINECOLOR_LABEL.SetEnabled(0);
btnLINECOLOR.SetEnabled(0);
btnACCEPT.SetEnabled(0);
btnCANCEL.SetEnabled(0);
GroupRemoveLayer(firstGroup, intersectVlayer);
# Delete the old station information and add the new station information
deleteStation();
addRecordToDatabase( addPointToVector(ctrPt) );
CloseVector(intersectV);
waitingForAddConfirmation = 0;
# Update the display
UpdateQuery();
ViewRedrawLayer(View, SDvectorLayer);
# Reset the polytool
tool.Managed = 0;
currentlyEditing = 0;
selectedElementNum = 0;
rGpMODE.SetEnabled(1);
btnCLOSE.SetEnabled(1);
promptString.SetLabel("Select point to Edit");
# Else the point needs to be loaded and allowed to be modified
} else {
promptString.SetLabel("Loading Data...");
btnACCEPT.SetEnabled(0);
btnCANCEL.SetEnabled(0);
overturnedToggle.SetEnabled(1);
aLINECOLOR_LABEL.SetEnabled(1);
btnLINECOLOR.SetEnabled(1);
# Get the original 3 points
local class POINT2D p1, p2, p3;
p1.x = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos1_x;
p1.y = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos1_y;
p2.x = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos2_x;
p2.y = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos2_y;
p3.x = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos3_x;
p3.y = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos3_y;
# Map them to the screen coordinates
local class TRANSPARM trans;
trans = ViewGetTransMapToView(View, rasterLayer.MapRegion.CoordRefSys , 0);
p1 = TransPoint2D(p1, trans);
p2 = TransPoint2D(p2, trans);
p3 = TransPoint2D(p3, trans);
trans = ViewGetTransViewToScreen(View, 0);
p1 = TransPoint2D(p1, trans);
p2 = TransPoint2D(p2, trans);
p3 = TransPoint2D(p3, trans);
# Reset the polytool with the 3 points
poly.Clear();
poly.AppendVertex(p1);
poly.AppendVertex(p2);
poly.AppendVertex(p3);
# Set the polytool for editing
currentlyEditing = 1;
tool.SetPolygon(poly);
tool.Managed = 1;
promptString.SetLabel("Modify the 3 points");
}
}
}
##################################################################################
#### Called when the user clicks on the Cancel button in the dialog. Actions are
#### based on the current Mode of the ToolScript
proc OnCancel( ) {
strikeFld.SetValue(0, 0);
dipFld.SetValue(0, 0);
dipDirFld.SetValue(0, 0);
# If in VIEW mode, unselect the point and reset the control dialog
if (Current_Mode == "VIEW") {
btnACCEPT.SetEnabled(0);
btnCANCEL.SetEnabled(0);
SDpoints.HighlightSingle(selectedElementNum, "Toggle");
selectedElementNum = 0;
rGpMODE.SetEnabled(1);
btnCLOSE.SetEnabled(1);
promptString.SetLabel("Select point to Toggle Data ON/OFF");
# If in ADD mode, remove the temporary outcrop trace pattern and reset the dialog
} else if (Current_Mode == "ADD") {
overturnedToggle.SetEnabled(0);
aLINECOLOR_LABEL.SetEnabled(0);
btnLINECOLOR.SetEnabled(0);
btnACCEPT.SetEnabled(0);
btnCANCEL.SetEnabled(0);
promptString.SetLabel("Clearing Data");
# Remove the computed temporary outcrop trace pattern and refresh the display
GroupRemoveLayer(firstGroup, intersectVlayer);
CloseVector(intersectV);
waitingForAddConfirmation = 0;
ViewRedrawLayer(View, SDvectorLayer);
ViewRedrawLayer(View, OPvectorLayer);
# Reset the tool
poly.Clear();
tool.SetPolygon(poly);
tool.Managed = 1;
rGpMODE.SetEnabled(1);
btnCLOSE.SetEnabled(1);
promptString.SetLabel("Select 3 new points");
# If in EDIT mode, cancel out of selecting a point, or remove the temporary data
} else if (Current_Mode == "EDIT") {
btnACCEPT.SetEnabled(0);
btnCANCEL.SetEnabled(0);
SDpoints.HighlightSingle(selectedElementNum, "Toggle");
tool.Managed = 0;
poly.Clear();
# Remove the temporary outcrop trace pattern data
if (currentlyEditing) {
GroupRemoveLayer(firstGroup, intersectVlayer);
CloseVector(intersectV);
}
currentlyEditing = 0;
waitingForAddConfirmation = 0;
rGpMODE.SetEnabled(1);
btnCLOSE.SetEnabled(1);
promptString.SetLabel("Select point to Edit");
# If in DELETE mode, deselect the point
} else if (Current_Mode == "DELETE") {
btnDELETE.SetEnabled(0);
btnCANCEL.SetEnabled(0);
SDpoints.HighlightSingle(selectedElementNum, "Toggle");
selectedElementNum = 0;
rGpMODE.SetEnabled(1);
btnCLOSE.SetEnabled(1);
promptString.SetLabel("Select point to Delete");
}
# Reset the selected element to nothing
selectedElementNum = 0;
}
##################################################################################
#### Called when the user clicks on the Cancel button in the dialog. Deletes the
#### selected point and all the records attached to that station.
proc OnDelete( ) {
promptString.SetLabel("Deleting Data...");
rGpMODE.SetEnabled(0);
btnCLOSE.SetEnabled(0);
btnCANCEL.SetEnabled(0);
btnDELETE.SetEnabled(0);
# Delete the station and attached records and update the display
deleteStation();
UpdateQuery();
ViewRedrawLayer(View, SDvectorLayer);
strikeFld.SetValue(0, 0);
dipFld.SetValue(0, 0);
dipDirFld.SetValue(0, 0);
rGpMODE.SetEnabled(1);
btnCLOSE.SetEnabled(1);
promptString.SetLabel("Select point to Delete");
selectedElementNum = 0;
}
##################################################################################
#### Called when the user clicks on the Close button to close the dialog
proc OnClose( ) {
tool.Managed = 0;
setDefaultWhenClose = true;
View.SetDefaultTool();
}
##################################################################################
#### Called when the user clicks on the "Set as Overturned" togglebutton. It
#### adjusts the value of the Strike angle by 180 degrees if it is set as overturned
proc OnOverturnedChanged( ) {
# If the Strike Angle is not 0
if (strikeFld.GetValueNum() != 0) {
# If the angle is less then 180 degrees and is being set to overturned
if ( (strikeAng < 180) && (overturnedToggle.GetValue()) ) {
strikeAng = strikeAng + 180;
# If the agnle is already overturned and being set as not overturned
} else if ( (strikeAng > 180) && (!overturnedToggle.GetValue()) ) {
strikeAng = strikeAng - 180;
}
# dlg.GetCtrlByID("id_STRIKEANGLE").SetValueNum(strikeAng);
strikeFld.SetValue(strikeAng,0);
}
}
##################################################################################
#### Called when the Right Button is pressed. This is an alternative method of
#### clicking "Accept" in certain situations
proc OnRightButtonPress( ) {
if (!tool.Managed) {
if (selectedElementNum > 0) {
if (Current_Mode == "VIEW") {
OnAccept();
} else if (Current_Mode == "DELETE") {
OnDelete();
}
}
}
}
##################################################################################
#### Called when the Left Button is pressed. This selects the closest point in the
#### StrikeDipVector and displays the information related to that point. Certain
#### modes use this selected point for their operations
proc OnLeftButtonPress( ) {
# Make sure the polytool is not currently active or it will conflict with it
if (!tool.Managed) {
# Get the point on the screen that was clicked
local class POINT2D screenPoint, layerPoint;
screenPoint.x = PointerX;
screenPoint.y = PointerY;
# Get the layer coordinates from the screen point and find the closest point
local class TRANSPARM transparm = ViewGetTransLayerToScreen(View, SDvectorLayer, 1);
layerPoint = TransPoint2D(screenPoint, transparm);
selectedElementNum = FindClosestPoint(StrikeDipVector, layerPoint.x, layerPoint.y, SDgeoref, 300);
# If greater then 0, it has found the closest point
if (selectedElementNum > 0) {
SDpoints.HighlightSingle(selectedElementNum, "Replace");
# Load the values for that point in the dialog
strikeFld.SetValue(StrikeDipVector.point[selectedElementNum].Bedding.StrikeAngle, 0);
dipFld.SetValue(StrikeDipVector.point[selectedElementNum].Bedding.DipAngle, 0);
dipDirFld.SetValue(StrikeDipVector.point[selectedElementNum].Bedding.DipDirection, 0);
overturnedToggle.SetEnabled(1);
overturnedToggle.SetValue(StrikeDipVector.point[selectedElementNum].Bedding.Overturned, 0);
overturnedToggle.SetEnabled(0);
# Find the color attached to that station and update the color button
local class COLOR tempColor;
local numeric selectedStation = StrikeDipVector.point[selectedElementNum].Bedding.Station;
local numeric numPoints = NumVectorPoints(OutcropPVector);
local numeric i = 1;
local numeric done = 0;
while (!done) {
# If it can find the station, load the color values
if (OutcropPVector.point[i].Outcrop.Station == selectedStation) {
tempColor.red = OutcropPVector.point[i].Outcrop.RedComp;
tempColor.green = OutcropPVector.point[i].Outcrop.GreenComp;
tempColor.blue = OutcropPVector.point[i].Outcrop.BlueComp;
done = 1;
# If not, load the color as Black and Notify user
} else if (i == numPoints) {
PopupMessage("Station Color Not Found");
tempColor.red = 0;
tempColor.green = 0;
tempColor.blue = 0;
done = 1;
}
i++;
}
btnLINECOLOR.SetColor(tempColor);
# Update the control dialog settings based on the current mode
if (Current_Mode == "VIEW") {
promptString.SetLabel("CONFIRM: Toggle Data for this point?");
rGpMODE.SetEnabled(0);
btnACCEPT.SetEnabled(1);
btnCANCEL.SetEnabled(1);
btnCLOSE.SetEnabled(0);
} else if (Current_Mode == "EDIT") {
promptString.SetLabel("CONFIRM: Edit this point?");
rGpMODE.SetEnabled(0);
btnACCEPT.SetEnabled(1);
btnCANCEL.SetEnabled(1);
btnCLOSE.SetEnabled(0);
} else if (Current_Mode == "DELETE") {
promptString.SetLabel("CONFIRM: Delete this point?");
rGpMODE.SetEnabled(0);
btnCANCEL.SetEnabled(1);
btnDELETE.SetEnabled(1);
btnCLOSE.SetEnabled(0);
}
}
}
}
###########################################################
### Procedure to check georeference of DEM; if it uses a geographic coordinate
### system, set up a coordinate transformation to the appropriate UTM zone
### to transform map coordinate used in determine strike/dip values.
### Called by onInitialize() and by cbLayer()
proc checkGeoref () {
local class RECT3D extentsObj;
local class POINT2D center;
transGeogCoords = 0;
demZscale = DEM.$Info.DataScale;
rvcGeorefDEM.OpenLastUsed(DEM);
demCRS = rvcGeorefDEM.GetCoordRefSys();
demCalibModel = rvcGeorefDEM.GetCalibModel();
rvcGeorefDEM.GetTransParm(transDEMobjToMap, 0, demCalibModel);
if (demCRS.IsProjected() <> 1) { # DEM CRS is not projected
if (demCRS.IsLocal() <> 1) { # DEM CRS also is not local, then must be geographic (lat/lon)
transGeogCoords = 1;
DEM.GetExtents(extentsObj);
center = extentsObj.center;
center = transDEMobjToMap.ConvertPoint2DFwd(center); # DEM center point in map coordinates, used to set UTM zone
# variables for setting up UTM CRS for transforming DEM map coordinates in tool
local class SR_COORDSYS utmCoordSys; # coordinate system
local class SR_COORDOPDEF projectDef; # projection parameters and method
local class SR_COORDOPMETHOD method;
local class SR_DATUM datum;
datum = demCRS.Datum; # get datum from the DEM raster
utmCoordSys.Assign("Projected2D_EN_m"); # assign projected coordinate system
method.Assign("1909"); # assign using numeric ID for Transverse Mercator projection
projectDef.Method = method; # assign method to projection
# assign projection parameters that are the same for all UTM zones
projectDef.SetParmValueNum("LatitudeOfNaturalOrigin", 0);
projectDef.SetParmValueNum("ScaleFactorAtNaturalOrigin", 0.9996);
projectDef.SetParmValueNum("FalseEasting", 500000);
# assign false northing based on whether north or south of equator
if ( center.y < 0) then
projectDef.SetParmValueNum("FalseNorthing", 10000);
else
projectDef.SetParmValueNum("FalseNorthing", 0);
# assign longitude of natural origin (central meridian of UTM zone) based on longitude
local numeric zoneNum, centMerid; # UTM zone and longitude of central meridian
zoneNum = ceil( (180 + center.x) / 6 );
centMerid = (zoneNum * 6) - 183;
projectDef.SetParmValueNum("LongitudeOfNaturalOrigin", centMerid);
# create defined UTM coordinate reference system
utmCRS.Create(utmCoordSys, datum, projectDef);
# setup coordinate transformation from DEM's geographic CRS to the defined UTM CRS.
transGeogUTM.InputCoordRefSys = demCRS;
transGeogUTM.OutputCoordRefSys = utmCRS;
}
}
}
##################################################################################
#### Called the first time the tool is activated. The dialog is created here, but
#### not yet displayed
proc OnInitialize ( ) {
# Is this a layout? Set callbacks in case things change
if (Layout) {
WidgetAddCallback(Layout.GroupSelectedCallback, cbGroup);
firstGroup = Layout.FirstGroup;
} else {
firstGroup = Group;
}
WidgetAddCallback(firstGroup.LayerSelectedCallback, cbLayer);
# Is the bottom layer a valid DEM raster?
if( checkLayer() ) {
# Initialize tool
tool = ViewCreatePolyLineTool(View);
ToolAddCallback(tool.ActivateCallback, cbToolApply);
checkGeoref();
# Create the destination vectors and their tables
createDestVectors();
addVectorsToDisplay();
initializeBeddingDatabaseTable();
initializeOutcropDatabaseTable();
View.Redraw();
} else {
string message$ = "Active Layer: "+activeLayer.Name+" is invalid.";
PopupMessage(message$);
}
# XML script that specifies the context of the dialog
string xml$='<?xml version="1.0"?>
<root>
<dialog id="id_StrikeDipDialog" Title="StrikeDip Toolscript Control Panel" buttons="">
<groupbox Name=" Current Mode: " ExtraBorder="5">
<radiogroup id="id_MODE" orientation="Horizontal" Default="ADD" OnSelection="OnModeChange();">
<item Value="VIEW" Name="View"/>
<item Value="ADD" Name="Add"/>
<item Value="EDIT" Name="Edit"/>
<item Value="DELETE" Name="Delete"/>
</radiogroup>
</groupbox>
<pane Orientation="Horizontal">
<label>Strike Azimuth: </label>
<editnumber Width="6" id="id_STRIKEANGLE" Precision="0" ReadOnly="true" HorizResize="Fixed"/>
</pane>
<pane Orientation="Horizontal">
<label>Dip Angle: </label>
<editnumber Width="6" id="id_DIPANGLE" Precision="0" ReadOnly="true" HorizResize="Fixed"/>
</pane>
<pane Orientation="Horizontal">
<label>Dip Azimuth: </label>
<editnumber Width="6" id="id_DIPDIRECTION" Precision="0" ReadOnly="true" HorizResize="Fixed"/>
</pane>
<togglebutton id="id_OVERTURNED" name="Set as Overturned" OnChanged="OnOverturnedChanged();" Selected="false" Enabled="false"/>
<pane Orientation="Horizontal">
<label id="id_LINECOLOR_LABEL">Set Outcrop Trace Color: </label>
<colorbutton id="id_LINECOLOR"/>
</pane>
<groupbox Name=" Apply Action: " Orientation="Horizontal" ExtraBorder="5">
<pushbutton id="id_ACCEPT" Name="Accept" OnPressed="OnAccept();"/>
<pushbutton id="id_CANCEL" Name="Cancel" OnPressed="OnCancel();"/>
<pushbutton id="id_DELETE" Name="Delete" OnPressed="OnDelete();"/>
<pushbutton id="id_CLOSE" Name="Close" OnPressed="OnClose();"/>
</groupbox>
<label id="id_PROMPT">Select Prompt</label>
</dialog>
</root>';
# Parse the document to make sure it is valid
local numeric err = dlgdoc.Parse(xml$);
if ( err < 0 ) {
PopupError( err ); # Popup an error dialog. "Details" button shows syntax errors.
Exit();
}
# Find the document
node = dlgdoc.GetElementByID("id_StrikeDipDialog");
if ( node == 0 ) {
PopupMessage("Could not find dialog node in XML document");
Exit();
}
# Set the XML Node
err = dlg.SetXMLNode(node);
if ( err < 0 ) {
PopupMessage("Could not set the XML Node");
Exit();
}
# Create a modeless dialog and assign some commonly changed elements to controls
dlg.CreateModeless();
rGpMODE = dlg.GetCtrlByID("id_MODE");
btnLINECOLOR = dlg.GetCtrlByID("id_LINECOLOR");
aLINECOLOR_LABEL = dlg.GetCtrlByID("id_LINECOLOR_LABEL");
overturnedToggle = dlg.GetCtrlByID("id_OVERTURNED");
promptString = dlg.GetCtrlByID("id_PROMPT");
strikeFld = dlg.GetCtrlByID("id_STRIKEANGLE");
dipFld = dlg.GetCtrlByID("id_DIPANGLE");
dipDirFld = dlg.GetCtrlByID("id_DIPDIRECTION");
btnACCEPT = dlg.GetCtrlByID("id_ACCEPT");
btnCANCEL = dlg.GetCtrlByID("id_CANCEL");
btnDELETE = dlg.GetCtrlByID("id_DELETE");
btnCLOSE = dlg.GetCtrlByID("id_CLOSE");
}
##################################################################################
#### Called when the tool is destroyed. It cleans up resources used.
proc OnDestroy ( ) {
tool.Managed = 0;
CloseDatabase(SDdb);
CloseDatabase(OPpointdb);
CloseDatabase(OPlinedb);
CloseVector(StrikeDipVector);
CloseVector(OutcropPVector);
}
##################################################################################
#### Called when tool is activated. It opens the dialog.
proc OnActivate ( ) {
dlg.Open();
OnModeChange();
setDefaultWhenClose = true;
}
##################################################################################
#### Called when tool is deactivated. Closes the dialog.
proc OnDeactivate ( ) {
dlg.Close(0);
setDefaultWhenClose = false;
}
Download