Syntax Highlighing:
comments, key words, predefined symbols, class members & methods, functions & classes
# REGSTATp.SML - Allows user to select an area of a non-composite raster
# by selecting and deselecting polygons from a vector layer over the top of it.
# The script then outputs the number of cells, number of non-null cells,
# minimum, maximum, mean, standard deviation, area, perimeter, centroid,
# and surface area over that region. The user may choose any distance or
# area units for the output to be displayed in.
#
# Requires TNTmips version 2012:78
# Version 26 August 2013
# Changed the function to create the pooint tool to ViewCreatePointGadget()
# to be consistent with internal class changes.
# The following symbols are predefined
# class VIEW View {use to access the view the tool script is attached to}
# class GROUP Group {use to access the group being viewed if the script is run from a group view}
# Variable declarations
class XmForm form, buttonRow, plusMinus;
class XmSeparator line1, line2;
class GRE_LAYER_VECTOR vectorLayer;
class GRE_LAYER rasterLayer;
class Vector targetVector;
class Raster targetRaster;
class XmDrawingArea da;
class GC gc;
class PointTool pointTool;
string vectorName$, rasterName$, regionMode$;
class PushButtonItem saveButton, closeButton;
class XmPushButton clearicon;
class ToggleButtonItem plusicon, minusicon;
class XmOptionMenu distMenu, areaMenu;
class REGION2D reg;
numeric min, max, mean, stdDev, count, cells, sum, sumsqr, surface, area, perimeter;
class POINT2D centroid;
numeric setDefaultWhenClose = false;
numeric distScale, areaScale, numPolys;
array numeric polyarray[1000];
class GRE_GROUP group;
# Checks layer to see if it is valid.
func checkLayer() {
local numeric valid = false;
# Get names layers if usable. If not output error messages.
if (group.ActiveLayer.TypeID == "") {
rasterName$ = "Group has no layers!";
}
else {
if (group.LastLayer.TypeID == "Vector") {
vectorLayer = group.LastLayer;
DispGetVectorFromLayer(targetVector, vectorLayer);
if (targetVector.$Info.NumPolys > 0) {
vectorName$ = vectorLayer.Name;
valid = true;
}
else vectorName$ = "No polygons!";
}
else vectorName$ = "Not a vector!";
if (group.FirstLayer.TypeID == "Raster") {
rasterLayer = group.FirstLayer;
DispGetRasterFromLayer(targetRaster, rasterLayer);
if (targetRaster.$Info.Type != "binary" and
targetRaster.$Info.Type != "4-bit unsigned" and
targetRaster.$Info.Type != "8-bit signed" and
targetRaster.$Info.Type != "8-bit unsigned" and
targetRaster.$Info.Type != "16-bit signed" and
targetRaster.$Info.Type != "16-bit unsigned" and
targetRaster.$Info.Type != "32-bit signed" and
targetRaster.$Info.Type != "32-bit unsigned" and
targetRaster.$Info.Type != "32-bit floating-point" and
targetRaster.$Info.Type != "64-bit signed" and
targetRaster.$Info.Type != "64-bit unsigned" and
targetRaster.$Info.Type != "64-bit floating-point") {
rasterName$ = "Type not supported!";
valid = false;
}
else
rasterName$ = rasterLayer.Name;
}
else {
rasterName$ = "Not a raster!";
valid = false;
}
}
return valid;
}
# Callback for drawing area expose. Draws text.
proc cbRedraw() {
local numeric larea, lperimeter, lsurface;
# Scale the values to specified units.
larea = areaScale * area;
lperimeter = distScale * perimeter;
lsurface = areaScale * surface;
if (gc == 0) return;
ActivateGC(gc);
# Clear the drawing area and redraw text.
SetColorName("gray75");
FillRect(0, 0, da.width, da.height);
SetColorName("black");
if (cells > 0) {
DrawInterfaceText(sprintf("Vector: %s\nRaster: %s\nCells: %d\nNull Cells: %d\nMinimum: %.2f\nMaximum: %.2f\nMean: %.2f\nStandard Deviation: %.2f\nArea: %.2f\nPerimeter: %.2f\nCentroid: %.2f, %.2f\nSurface Area: %.2f",vectorName$, rasterName$, count, cells - count, min, max, mean, stdDev, larea, lperimeter, centroid.x, centroid.y, lsurface), 0, 10);
}
else DrawInterfaceText(sprintf("Vector: %s\nRaster: %s\nCells:\nNull Cells:\nMinimum:\nMaximum:\nMean:\nStandard Deviation:\nArea:\nPerimeter:\nCentroid:\nSurface Area:",vectorName$, rasterName$), 0, 10);
}
# Callback for when the active layer changes.
proc cbLayer() {
checkLayer();
cbRedraw();
}
# Callback for when the active group changes.
proc cbGroup() {
group = Layout.ActiveGroup;
WidgetAddCallback(group.LayerSelectedCallback, cbLayer);
cbLayer();
}
# Callback for when user clicks the right mouse button on the polygon tool
# (or clicks apply).
proc cbToolApply(class pointTool tool) {
# If the selected layer is not valid, don't do anything.
if (checkLayer()) {
# Set local variables
local numeric sum, sumsqr, xscale, yscale, zscale;
local numeric current, right, down, downright, elemNum;
local class POINT2D point;
local class StatusHandle status;
local class StatusContext context;
cells = 0; min = 0; max = 0; mean = 0; stdDev = 0; sum = 0; sumsqr = 0; count = 0; surface = 0; area = 0; perimeter = 0;
centroid.x = 0; centroid.y = 0; current = 0; right = 0; down = 0; downright = 0;
xscale = ColScale(targetRaster);
yscale = LinScale(targetRaster);
zscale = group.FirstLayer.zscale;
# Check point.
point.x = tool.Point.x;
point.y = tool.Point.y;
point = TransPoint2D(point, ViewGetTransViewToScreen(View, 1));
point = TransPoint2D(point, ViewGetTransMapToView(View, vectorLayer.Projection, 1));
numeric elementNum = FindClosestPoly(targetVector, point.x, point.y, GetLastUsedGeorefObject(targetVector));
if (elementNum > 0) {
if (regionMode$ == "plus") {
# Add to polygon list.
numPolys += 1;
polyarray[numPolys] = elementNum;
vectorLayer.Poly.HighlightSingle(elementNum, "Add");
}
else {
if (numPolys > 1) {
# Remove from polygon list.
local numeric i;
for i = 1 to numPolys
if (polyarray[i] == elementNum) {
polyarray[i] = polyarray[numPolys];
numPolys -= 1;
vectorLayer.Poly.HighlightSingle(elementNum, "Subtract");
break;
}
}
else if (numPolys == 1 and polyarray[1] == elementNum) {
vectorLayer.Poly.HighlightSingle(elementNum, "Subtract");
numPolys = 0;
}
}
# Get region and transform it to the appropriate coordinate system.
if (numPolys > 0) {
reg = ConvertVectorPolysToRegion(targetVector, GetLastUsedGeorefObject(targetVector), polyarray, numPolys);
}
else ClearRegion(reg);
# Create a status bar so user knows the script is doing something.
status = StatusDialogCreate(form);
context = StatusContextCreate(status);
StatusSetMessage(context, "Computing values...");
# Loop on region to calculate statistics.
local numeric lin, col;
foreach targetRaster[lin, col] in reg {
if (!IsNull(targetRaster)) {
if (count == 0) {
max = targetRaster;
min = targetRaster;
}
else if (targetRaster > max) {
max = targetRaster;
}
else if (targetRaster < min) {
min = targetRaster;
}
sum += targetRaster;
sumsqr += sqr(targetRaster);
count += 1;
}
# Surface area is estimated by by adding the areas of triangles created
# by the current cell, the cell below it, the cell to the right of it,
# and the cell to the lower right of it.
# Estimate over null cells by keeping the last known value.
if (!IsNull(targetRaster))
current = targetRaster;
if (!IsNull(targetRaster[lin,col+1]))
right = targetRaster[lin,col+1];
if (!IsNull(targetRaster[lin+1,col]))
down = targetRaster[lin+1,col];
if (!IsNull(targetRaster[lin+1,col+1]))
downright = targetRaster[lin+1,col+1];
surface += .5*sqrt(sqr(yscale*current*zscale-yscale*right*zscale)
+sqr(xscale*current*zscale-xscale*down*zscale)
+sqr(xscale*yscale));
surface += .5*sqrt(sqr(yscale*downright*zscale-yscale*right*zscale)
+sqr(xscale*downright*zscale-xscale*down*zscale)
+sqr(xscale*yscale));
cells += 1;
}
CloseRaster(targetRaster);
# Avoid division by zero when calculating mean and standard deviation.
if (count > 1) {
mean = sum / count;
stdDev = sqrt((sumsqr - sqr(sum) / count) / (count - 1));
area = reg.GetArea();
perimeter = reg.GetPerimeter();
centroid = reg.GetCentroid();
}
cbRedraw();
# Destroy the status bar
StatusContextDestroy(context);
StatusDialogDestroy(status);
}
}
}
# Called when a unit type is selected.
proc cbDistUnits(class Widget parent, unit$) {
distScale = GetUnitConvDist("meters", unit$);
cbRedraw();
}
# Called when a unit type is selected.
proc cbAreaUnits(class Widget parent, unit$) {
areaScale = GetUnitConvArea("square meters", unit$);
cbRedraw();
}
# Called when the save button is pressed. Saves statistics to text file.
proc cbSave() {
local class FILE outfile;
local numeric larea, lperimeter, lsurface;
# Scale the values to specified units.
larea = areaScale * area;
lperimeter = distScale * perimeter;
lsurface = areaScale * surface;
# Open output file, write to it, and close it
outfile = GetOutputTextFile(_context.ScriptDir, "Choose file...", "txt");
fprintf(outfile, "Raster: %s\nCells: %d\nNull Cells: %d\nMinimum: %.2f\nMaximum: %.2f\nMean: %.2f\nStandard Deviation: %.2f\nArea: %.2f\nPerimeter: %.2f\nCentroid: %.2f, %.2f\nSurface Area: %.2f",
rasterName$, count, cells - count, min, max, mean, stdDev, larea, lperimeter, centroid.x, centroid.y, lsurface);
fprintf(outfile, "\nDistance Units: %s\nArea Units: %s\n\n", distMenu.value, areaMenu.value);
fclose(outfile);
}
# Called when the user presses the clear region button.
proc cbClearRegion() {
vectorLayer.UnHighlightAllElements();
numPolys = 0;
}
# Called when the close button is pressed. Closes the dialogs.
proc cbClose() {
pointTool.Managed = 0;
pointTool.HasPosition = 0;
DialogClose(form);
if (setDefaultWhenClose) {
setDefaultWhenClose = false;
View.SetDefaultTool();
}
}
# Called when the user selects plus region.
proc cbPlusRegion() {
regionMode$ = "plus";
}
# Called when the user selects minus region.
proc cbMinusRegion() {
regionMode$ = "minus";
}
# Called the first time the tool is activated.
# If the tool implements a dialog it should be created (but not displayed) here.
func OnInitialize () {
if (Layout) {
WidgetAddCallback(Layout.GroupSelectedCallback, cbGroup);
group = Layout.ActiveGroup;
}
else
group = Group;
WidgetAddCallback(group.LayerSelectedCallback, cbLayer);
# Set up tool dialog.
form = CreateFormDialog("Region Statistics");
form.marginHeight = 2;
form.marginWidth = 2;
WidgetAddCallback(form.Shell.PopdownCallback, cbClose);
# add region mode (+,-) region icons
# set the icons to initial state - plus region
plusicon = CreateToggleButtonItem("Add to region", cbPlusRegion);
plusicon.IconClass = "standard";
plusicon.IconName = "add_select";
plusicon.OneOfMany = true;
minusicon = CreateToggleButtonItem("Subtract from region", cbMinusRegion);
minusicon.IconClass = "standard";
minusicon.IconName = "minus_select";
minusicon.OneOfMany = true;
plusMinus = CreateForm(form);
CreateIconButtonRow(plusMinus, plusicon, minusicon);
plusMinus.TopWidget = form;
plusMinus.LeftWidget = form;
clearicon = CreateIconPushButton(form, "standard", "delete", "Clear regions");
clearicon.TopWidget = form;
clearicon.leftWidget = plusMinus;
WidgetAddCallback(clearicon.ActivateCallback, cbClearRegion);
# Drawing area for displaying text.
da = CreateDrawingArea(form, 185, 301);
da.topWidget = plusMinus;
da.leftWidget = form;
da.rightWidget = form;
WidgetAddCallback(da.ExposeCallback, cbRedraw);
# Separator between output and unit tools.
line1 = CreateHorizontalSeparator(form);
line1.topWidget = da;
line1.leftWidget = form;
line1.rightWidget = form;
line1.bottomOffset = 2;
# Menus for selecting units.
distMenu = CreateUnitOptionMenu(form, "distance_units_c", cbDistUnits, 2, 0);
distMenu.topWidget = line1;
distMenu.leftWidget = form;
areaMenu = CreateUnitOptionMenu(form, "area_units_c", cbAreaUnits, 1, 0);
areaMenu.topWidget = distMenu;
areaMenu.leftWidget = form;
# Separator between unit tools and buttons.
line2 = CreateHorizontalSeparator(form);
line2.topWidget = areaMenu;
line2.leftWidget = form;
line2.rightWidget = form;
line2.topOffset = 2;
# Allows user to save statistics to text file.
saveButton = CreatePushButtonItem("Save As...", cbSave);
# Closes the dialog.
closeButton = CreatePushButtonItem("Close", cbClose);
# Row of buttons.
buttonRow = CreateButtonRow(form, saveButton, closeButton);
buttonRow.topWidget = line2;
buttonRow.leftWidget = form;
buttonRow.rightWidget = form;
# Add point tool.
pointTool = ViewCreatePointGadget(View);
ToolAddCallback(pointTool.ActivateCallback, cbToolApply);
} # end of OnInitialize
# Called when tool is to be destroyed, will not be called if tool was never activated.
func OnDestroy () {
pointTool.Managed = 0;
DestroyGC(gc);
DestroyWidget(form);
} # end of OnDestroy
# Called when tool is activated.
func OnActivate () {
pointTool.Managed = 1;
pointTool.HasPosition = 0;
plusicon.Selected = 1;
if (checkLayer())
cbClearRegion();
DialogOpen(form);
if (gc == 0)
gc = CreateGCForDrawingArea(da);
cells = 0; min = 0; max = 0; mean = 0; stdDev = 0; sum = 0; sumsqr = 0; count = 0; surface = 0; area = 0; perimeter = 0;
centroid.x = 0; centroid.y = 0;
distScale = 1;
areaScale = 1;
distMenu.value = "meters";
areaMenu.value = "square meters";
cbRedraw();
setDefaultWhenClose = true;
} # end of OnActivate
# Called when tool is deactivated (usually when switching to another tool).
func OnDeactivate () {
setDefaultWhenClose = false;
cbClose();
} # end of OnDeactivate