#
# Biomass.sml (SML object within the biomass2.sml Applidat)
#
# This is how to declare a global variable so that
# they can be used in functions before it's assigned.
# Without the explicit declaration, the variable will
# be assumed to be numeric or a Raster (depending on
# the case of the first letter).
#
numeric sample_height, sample_width, sample_space, sample_text_x, num_ranges;
sample_height = 15; # Height of sample boxes
sample_width = 40; # Width of sample boxes
sample_space = 5; # Space between boxes
sample_text_x = 37; # Where to start the right justified text.
num_ranges = 10; # Number of ranges to break histogram into.
class ColorMap cmap;
class Color black, white, c0;
class XmDrawingArea da;
class XmForm form; # The main form dialog
class GC gc; # The context for drawing into the view
class RASTER Composite, BiomassRaster, FilteredBiomassRaster;
class VECTOR BiomassVector;
class GRE_VIEW view;
class GRE_GROUP group;
class GRE_LAYER_RASTER layer, overlay; # Our base layer
class GRE_LAYER_VECTOR voverlay;
array numeric histo[256];
numeric cellsize, unitconv;
class GUI_GADGET_POLYLINE tool;
class GUI_CTRL_TOGGLEBUTTON plusbutton;
class GUI_CTRL_TOGGLEBUTTON minusbutton;
class XmToggleButton drawhisto;
class PushButtonItem buttonConvertToVector, buttonFiltVector, buttonFiltRaster;
class PushButtonItem buttonClearOverlay, buttonExit;
class XmForm button_row;
class XmFrame frame;
numeric height, width;
class REGION2D MyRgn;
class ToolTip tip, tip2, tip3;
string str$ = "";
numeric firstpass, minusregion;
firstpass = 1;
minusregion = 0;
#unitconv = GetUnitConvArea("square meters", "acres");
# Note, GetUnitConvArea() fails if running on a system with localized
# resources. So just hard code the conversion
unitconv = 1.0 / 4046.873;
# If the script is in an RVC file, assume all the data
# is in that file. If it's not, assume I'm still working
# on it and get the data from my data directory
if (_context.ScriptIsRVCObject) {
string datafile$ = _context.Filename;
}
else {
datafile$ = "e:/data/toolbar/smldemo.rvc";
}
string outputfile$ = FileNameGetPath(datafile$) + "/biomass.rvc";
CreateProjectFile(outputfile$, "Biomass");
################################################################
#
# This is the callback function which will be called
# when the exit button is pressed.
#
proc cbQuit(class widget widget) {
DialogClose(form);
}
################################################################
#
# This is the callback function which will be called
# when the dialog goes away. This can happen either because
# we called DialogClose() from our Quit callback or because
# the user closed the dialog using the [X] icon on the title bar.
#
proc cbPopdown(class widget widget) {
ViewDestroy(view);
GroupDestroy(group);
DestroyWidget(form);
Exit(); # Nothing else from here on!
}
################################################################
#
# Draw the bar graph
#
proc drawGraph(class XmDrawingArea da) {
local numeric maxval, total, i, scale, w, y;
# Clear the drawing area
if (!gc) return;
ActivateGC(gc);
SetColorRGB(0,0,0);
FillRect(0, 0, da.width, da.height);
maxval = 0;
total = 0;
for i = 1 to num_ranges {
maxval = SetMax(maxval, histo[i]);
total += histo[i];
}
if (total) { # Don't draw bars if nothing selected
w = da.width - 2 * sample_space;
scale = w / maxval;
if (drawhisto.set) {
for i = 1 to num_ranges {
y = (num_ranges - i + 1) * (sample_height + sample_space) + sample_space;
SetColor(ColorMapGetColor(cmap, i));
w = scale * histo[i];
if (w < 1 && histo[i] != 0) w = 1;
FillRect(sample_space, y, w, sample_height);
}
}
else {
DrawTextSetColors(white, black);
y = sample_space;
gc.TextStyle.JustifyRight = 1;
DrawTextSetColors(black, white);
for i = 1 to num_ranges {
y = (num_ranges - i) * (sample_height + sample_space) + sample_space;
SetColor(ColorMapGetColor(cmap, i));
FillRect(sample_space, y, sample_width, sample_height);
SetColorName("white");
DrawRect(sample_space, y, sample_width, sample_height);
str$ = sprintf("%.2f", histo[i] * cellsize * unitconv);
DrawTextSimple(str$, sample_text_x, y + sample_height - 4);
}
y = (num_ranges) * (sample_height + sample_space) + sample_space;
SetColorName("white");
MoveTo(2, y);
DrawTo(da.width - 2, y);
str$ = sprintf("%.2f", total * cellsize * unitconv);
DrawTextSetColors(white, black);
y += (sample_height - 3);
DrawTextSimple(str$, sample_text_x, y);
DrawTextSimple("TOTAL", sample_text_x, y + 11);
DrawTextSimple("ACRES", sample_text_x, y + 22);
}
}
}
################################################################
#
# Callback for when user clicks the right mouse button
# on the polygon tool (or hits apply)
#
proc cbToolApply(class GUI_GADGET_POLYLINE tool) {
local numeric i, val, scale, minval, maxval, ir, red, total, offset, lo, hi;
local class REGION2D reg, reg2;
class POLYLINE poly = tool.GetPolygon();
poly.SetClosed(1);
reg.UnionPolyLine(poly);
if (firstpass) {
MyRgn = CopyRegion(reg);
firstpass = 0;
}
else {
if (minusregion) {
reg2 = RegionSubtract(MyRgn, reg);
MyRgn = CopyRegion(reg2);
reg = CopyRegion(reg2);
}
else {
# reg = RegionOR(MyRgn, tool.Region);
reg2 = RegionOR(MyRgn, reg);
MyRgn = CopyRegion(reg2);
reg = CopyRegion(reg2);
}
}
PopupMessage(NumToStr(reg.Extents.x1) + " " + NumToStr(reg.Extents.x2));
reg = RegionTrans(reg,ViewGetTransLayerToScreen(view, layer, 1));
StatusSetDefaultHandle(view.StatusHandle);
ViewSetMessage(view, "Filling raster with 0");
BiomassRaster = 0;
if (1) {
ViewSetMessage(view, "Computing scale");
for i=0 to 200 histo[i] = 0;
for each Composite in reg {
ir = Composite.red;
red = Composite.green;
val = ((ir - red) / (ir + red) + 1) * 100;
histo[val] = histo[val] + 1;
}
total = 0;
for i=0 to 200 total += histo[i];
lo = total * .02;
hi = total * .98;
minval = -100;
maxval = -100;
total = 0;
for i = 0 to 200 {
if (total > lo && minval < -1) {
minval = (i / 100) - 1;
}
total += histo[i];
if (total > hi && maxval < -1) {
maxval = (i / 100) - 1;
}
}
scale = num_ranges / (maxval - minval);
offset = -minval;
}
else {
scale = num_ranges / 2.0;
offset = 1.0;
}
PopupMessage(NumToStr(reg.Extents.x1) + " " + NumToStr(reg.Extents.x2));
ViewSetMessage(view, "Processing region");
for i = 1 to num_ranges histo[i] = 0;
for each Composite in reg {
ir = Composite.red;
red = Composite.green;
val = (ir - red) / (ir + red);
val = Bound((val + offset) * scale + 1, 1, num_ranges);
histo[val] = histo[val] + 1;
BiomassRaster = val;
}
CloseRaster(BiomassRaster);
LayerShow(overlay, view);
ViewRedrawLayer(view, overlay);
ViewRedraw(view);
drawGraph(da);
# Do something with the region here
tool.HasPosition = 0; # Clears the current position
buttonConvertToVector.Disabled = 0;
buttonFiltRaster.Disabled = 0;
}
################################################################
#
# Callback for the "Filter Biomass" button
#
proc cbFiltRast () {
local class REGION2D reg;
reg = RegionTrans(MyRgn, ViewGetTransLayerToView(view, layer, 1));
FilteredBiomassRaster = 0;
ViewSetMessage(view, "Median filter");
StatusSetDefaultHandle(view.StatusHandle);
foreach BiomassRaster in reg {
FilteredBiomassRaster = FocalMajority(BiomassRaster, 5, 5);
}
foreach BiomassRaster in reg {
#BiomassRaster = FocalMajority(FilteredBiomassRaster, 3, 3);
BiomassRaster = (FilteredBiomassRaster);
}
CloseRaster(BiomassRaster);
ViewRedraw(view);
drawGraph(da);
}
################################################################
#
# Callback for the "Convert To Vector" button
#
proc cbConvToVector () {
local numeric i, numPolys, numToDelete, minsize, inode;
minsize = IniReadNumber("biomass", "minsize", 0);
minsize = PopupNum("Enter minimum polygon size in acres:", minsize, 0, 640, 1);
IniWriteNumber("biomass", "minsize", minsize);
minsize /= unitconv; # Convert to sq meters
if (VectorExists(outputfile$, "Vector")) {
inode = ObjectNumber(outputfile$, "Vector", "VECTOR");
DeleteObject(outputfile$, inode);
}
CreateVector(BiomassVector, outputfile$, "Vector", "Vector boundaries");
FilteredBiomassRaster = 0;
ViewSetMessage(view, "Median filter");
foreach BiomassRaster in MyRgn {
FilteredBiomassRaster = FocalMedian(BiomassRaster, 5, 5);
}
StatusSetDefaultHandle(view.StatusHandle);
ViewSetMessage(view, "Converting Raster to Vector");
BiomassVector = RasterToVectorBound(BiomassRaster);
ViewSetMessage(view, "Generating std stats table");
VectorToolkitInit(BiomassVector);
VectorUpdateStdAttributes(BiomassVector); # Compute std stats table
numToDelete = 0;
numPolys = NumVectorPolys(BiomassVector);
array numeric deleteList[numPolys];
ViewSetMessage(view, "Looking for polygons to delete");
for i=1 to numPolys {
if (BiomassVector.poly[i].POLYSTATS.Area < minsize) {
numToDelete += 1;
deleteList[numToDelete] = i;
}
}
if (numToDelete) {
ViewSetMessage(view, "Deleting tiny polygons");
VectorDeletePolys(BiomassVector, deleteList, numToDelete);
}
ViewSetMessage(view, "Closing vector");
CloseVector(BiomassVector);
if (!voverlay) {
voverlay = GroupQuickAddVectorVar(group, BiomassVector);
voverlay.Poly.NormalStyle.BorderColor.name = "black";
voverlay.Line.NormalStyle.Color.name = "black";
voverlay.Line.NormalStyle.width = 1;
voverlay.Line.NormalStyle.MapScale = 2000;
voverlay.Line.NormalStyle.ScaleToMap = 1;
}
LayerShow(voverlay, view);
ViewRedraw(view);
drawGraph(da);
}
proc cbHideOverlay(class widget widget) {
local numeric i;
tool.HasPosition = 0; # Clears the current position
buttonConvertToVector.Disabled = 1;
firstpass = 1;
buttonFiltRaster.Disabled = 1;
LayerHide(overlay, view);
if (voverlay) {
LayerHide(voverlay, view);
}
for i = 0 to 10 histo[i] = 0;
ViewRedraw(view);
drawGraph(da);
}
################################################################
#
# A generic callback to redraw the graph. Will be called
# after every redraw because the colormap may have changed
# out from under us.
#
proc cbRedrawViewBegin(class GRE_VIEW view) {
tip.Disabled = 1;
tip2.Disabled = 1;
tip3.Disabled = 1;
}
proc cbRedrawViewEnd(class GRE_VIEW view) {
# Force all the tips timers to restart
tip.Disabled = 1;
tip2.Disabled = 1;
tip3.Disabled = 1;
tip.Disabled = 0;
tip2.Disabled = 0;
tip3.Disabled = 0;
drawGraph(da);
}
proc cbRedrawGraph(class widget button) {
drawGraph(da);
}
proc cbMinus(class GUI_CTRL_TOGGLEBUTTON this) {
if (minusbutton.GetValue()) plusbutton.SetValue(0, false);
else plusbutton.SetValue(1, false);
minusregion = minusbutton.GetValue();
}
proc cbPlus(class GUI_CTRL_TOGGLEBUTTON this) {
if (plusbutton.GetValue()) minusbutton.SetValue(0, false);
else minusbutton.SetValue(1, false);
minusregion = minusbutton.GetValue();
}
#####################################################
#
# MAIN PROGRAM
#
# clear()
# Get an input raster. We will use it's size to determine
# how big to make our window
OpenRaster(Composite, datafile$, "CIR_Raster");
#GetInputRaster(Composite);
#-----------------------------------------------------
#
# Create the Biomass Raster that we'll use for the overlay
#
#CreateTempRaster(BiomassRaster, NumLins(Composite), NumCols(Composite), "8-bit unsigned", 1);
string lokfile$ = FileNameGetPath(outputfile$) +"/"+ FileNameGetName(outputfile$) + ".lok";
#DeleteFile(lokfile$);
DeleteFile(outputfile$);
CreateRaster(BiomassRaster, outputfile$, "biomass", "biomass raster", NumLins(Composite), NumCols(Composite), "8-bit unsigned");
CopySubobjects(Composite, BiomassRaster, "georef");
cellsize = LinScale(Composite) * ColScale(Composite);
#ini = IniOpenObject(datafile$, "Prefs");
IniWriteString("BiomassRaster", "Filename", BiomassRaster.$Info.Filename);
IniWriteString("BiomassRaster", "ObjectName", BiomassRaster.$Info.Name);
#IniClose(ini);
c0 = ColorMapGetColor(cmap, 0);
c0.Name = "black";
c0.transp = 75; # Make it 75% transparent
ColorMapSetColor(cmap, 0, c0);
ColorMapSetColorHIS(cmap, 10, 230, 50, 100); # Green
ColorMapSetColorHIS(cmap, 9, 180, 70, 100); # Yellow
ColorMapSetColorHIS(cmap, 8, 150, 70, 100); # Lt Orange
ColorMapSetColorHIS(cmap, 7, 144, 50, 100); # Dk Orange
ColorMapSetColorHIS(cmap, 6, 120, 50, 100); # Red
ColorMapSetColorHIS(cmap, 5, 72, 50, 100); # Magenta
ColorMapSetColorHIS(cmap, 4, 48, 50, 100); # Purple
ColorMapSetColorHIS(cmap, 3, 10, 50, 100); # Blue
ColorMapSetColorHIS(cmap, 2, 342, 60, 100); # Lt Blue
ColorMapSetColorHIS(cmap, 1, 320, 70, 100); # Cyan
#for i = 0 to 10 {
# c0 = ColorMapGetColor(cmap, i);
# c0.transp = 75; # Make it 75% transparent
# ColorMapSetColor(cmap, i, c0);
# }
ColorMapWriteToRastVar(BiomassRaster, cmap, "ColorMap");
SetNull(BiomassRaster, 0);
CreateHistogram(BiomassRaster);
CreateTempRaster(FilteredBiomassRaster, NumLins(Composite), NumCols(Composite), "8-bit unsigned");
CopySubobjects(Composite, FilteredBiomassRaster, "georef");
#-----------------------------------------------------
#
# Create a temp vector for the biomass->vector conversion
#
#CreateTempVector(BiomassVector);
#CreateVector(BiomassVector, outputfile$, "Vector", "Vector boundaries");
# We don't create it until we actually need it now
#-----------------------------------------------------
# Create a dialog (aka "window").
#
# This function returns a form that can be used as a
# parent for other widgets
#
# If we had another dialog which we wanted to use as
# a "parent" to this one, we would pass it as the
# 2nd parameter to CreateModalFormDialog().
form = CreateFormDialog("Biomass Mapping");
WidgetAddCallback(form.Shell.PopdownCallback, cbPopdown);
#-----------------------------------------------------
#
# Create a push button on the form. Set its attachments
# so that sets in the lower left of the form. To attach
# an edge to the form rather than another widget, just
# set it to the form. In Motif, this would be done by
# setting the attachment to XmATTACH_FORM, which I plan
# to implement here too. Setting the RightPosition to 50
# ties the right edge of the button to a spot 50% of the
# width of the form, (assuming form.FractionBase still
# has the default value of 100) no matter how wide the
# form gets
#
buttonClearOverlay = CreatePushButtonItem("Clear Overlay", cbHideOverlay);
buttonFiltRaster = CreatePushButtonItem("Filter Biomass", cbFiltRast);
buttonConvertToVector = CreatePushButtonItem("Convert To Vector", cbConvToVector);
#buttonFiltVector = CreatePushButtonItem("Filter Vector", cbFilterVector);
buttonExit = CreatePushButtonItem("Exit", cbQuit);
buttonConvertToVector.Disabled = 1;
buttonFiltRaster.Disabled = 1;
#buttonFiltVector.Disabled = 1;
button_row = CreateButtonRow(form, buttonClearOverlay,
buttonFiltRaster,
buttonConvertToVector,
buttonExit
);
button_row.BottomWidget = form;
drawhisto = CreateIconToggleButton(form, "RVCobjects", "histo", "Show as histogram");
drawhisto.LeftWidget = form;
drawhisto.TopWidget = form;
WidgetAddCallback(drawhisto.ValueChangedCallback, cbRedrawGraph);
frame = CreateFrame(form);
frame.TopWidget = drawhisto;
frame.BottomWidget = button_row;
frame.LeftWidget = form;
da = CreateDrawingArea(frame, 180, sample_width + 2 * sample_space + 2);
WidgetAddCallback(da.ExposeCallback, cbRedrawGraph);
#-----------------------------------------------------
#
# Create a 2D group to display the raster in and quick-add it.
#
group = DispCreate2DGroup();
#-----------------------------------------------------
#
# Create a view for this group in our form. Default its size to
# the size of the raster. For a real program, you'd want to
# limit this to something reasonable.
#
height = NumLins(Composite);
width = NumCols(Composite);
while (height > 700 or width > 1000) {
height /= 2;
width /= 2;
}
view = GroupCreateView(group, form, "view", height, width);
view.DisableRedraw = 1; #Prevent open from doing a redraw because
DispAddCallback(view.DrawBeginCallback, cbRedrawViewBegin);
DispAddCallback(view.DrawEndCallback, cbRedrawViewEnd);
DispAddCallback(view.RestoreCallback, cbRedrawViewEnd);
#-----------------------------------------------------
#
# Make widget attachments for the view. Connect it to the form
# on all edges except the bottom. There, attach it to one of the
# buttons.
# Some notes about making attachment:
# 1) If you attach widget A to Widget B, don't attach B to A
# or anything else which would cause a circular reference.
# 2) The widget created and attached last is what will resize
# with the window. If I had created the view first and then
# the buttons, resizing the window would make the buttons
# bigger.
form.TopWidget = form;
form.LeftWidget = frame;
form.RightWidget = form;
form.BottomWidget = button_row;
view.ShowDataTips = 0;
#-----------------------------------------------------
#
# Create a region tool on the view
#
tool = ViewCreatePolyLineTool(view, "My Tool", "draw_polygon", "polyline_edit");
ToolAddCallback(tool.ActivateCallback, cbToolApply);
tool.DialogPosition = "RightCenter";
#-----------------------------------------------------
#
# Add the standard graphic tools to the view
#
ViewAddStandardTools(view);
#-----------------------------------------------------
#
# Have to call this function after adding any tools to
# create icon buttons for the tools.
#
ViewAddToolIcons(view);
#-----------------------------------------------------
#
# add region mode (+,-) region icons
#
#plusicon = CreateToggleButtonItem("Add to region", cbPlusRegion);
#minusicon = CreateToggleButtonItem("Subtract from region", cbMinusRegion);
# set the icons to initial state - plus region
#plusicon.IconClass = "standard";
#plusicon.IconName = "add_select";
#plusicon.Selected = true;
#plusicon.OneOfMany = true;
#minusicon.IconClass = "standard";
#minusicon.IconName = "minus_select";
#minusicon.OneOfMany = true;
plusbutton.CreateIcon(view.ToolBarPane, "CONTROL_ADD_CYAN", "Add to region");
plusbutton.SetOnPressed(cbPlus);
plusbutton.SetValue(1, false);
minusbutton.CreateIcon(view.ToolBarPane, "CONTROL_SUBTRACT_CYAN", "Subtract from region");
minusbutton.SetOnPressed(cbMinus);
#-----------------------------------------------------
#
# Create the layers
#
layer = GroupQuickAddRasterVar(group, Composite);
layer.NullCellsTransparent = 1;
overlay = GroupQuickAddRasterVar(group, BiomassRaster);
overlay.NullCellsTransparent = 0;
overlay.UsesTransparency = 1;
overlay.DoBlendMask = 1;
LayerHide(overlay, view);
#-----------------------------------------------------
#
# Open the dialog Need to do this before creating the GC
# because until then the drawing area has no window.
# Would like a way around this. Could my CreateGC
# realize it first if necessary?
#
DialogOpen(form); # we're just going to do a resize right away
DialogFullScreen(form);
view.DisableRedraw = 0; # and that will make us redraw again.
ViewRedrawIfNeeded(view);
#ViewActivateTool(view, tool);
tool.Managed = 1;
#
# Create a GC (Graphics Context) for the drawing area
# You can have multiple gc's per drawing area. Drawing
# functions will always use the active GC.
#
gc = CreateGCForDrawingArea(da);
gc.TextStyle.JustifyRight = 1;
# Clear the drawing area
ActivateGC(gc);
SetColorRGB(0,0,0);
black.name = "black";
white.name = "white";
FillRect(0, 0, da.width, da.height);
DrawTextSetFont("stork.of");
DrawTextSetColors(black, white);
DrawTextSetHeightPixels(9);
# Add the layers after we've had a chance to notice the resize
# Don't do the following. Need the tool window to show.
#DialogToTop(form); # Pull our window to the top
#-----------------------------------------------------
#
# Add the data tips.
#
ViewSetMessage(view, 'Draw around a field and press "Apply"');
tip = CreateToolTip(view.DrawingArea,
'(1) Press the left mouse button to draw around the
field whose canopy biomass is of interest.
Click on "apply" to close your field boundary'
);
tip.Disabled = 1;
tip.Delay = 5000; # pop up after 3000 ms
tip.priority = -3; # Less than DataTips
tip2 = CreateToolTip(view.DrawingArea,
'(2) Press the left mouse button to draw around your field.
If you want to add another field or add area to a field
just draw it now. If you want to omit any area within a
field use left mouse button to select the "-" (minus)
icon in the upper right of the menu toolbar. You will
then be in the subtract area mode. Next simply draw
around the area to be deleted.'
);
tip2.Disabled = 1;
tip2.Delay = tip.Delay + 8000; # pop up after another 6000 ms
tip2.priority = -2; # Less than DataTips
tip3 = CreateToolTip(view.DrawingArea,
'(3) To redraw, zoom, pan... use the button icons on the left
of the menu toolbar. For more help use the "Instruction"
icon on the product toolbar.'
);
tip3.Disabled = 1;
tip3.Delay = tip2.Delay + 18000; # pop up after another 10 sec
tip3.priority = -1; # Less than DataTips
tip.Disabled = 0;
tip2.Disabled = 0;
tip3.Disabled = 0;
#-----------------------------------------------------
#
# New and improved method!
# Rather than having the script wait for the dialog
# to close, setup an OnExit callback. If the script
# is exited for any reason (including clicking the
# exit icon on the tool bar), this function will be
# called before we're terminated.
# Then just call WaitForExit(). This will put the
# script into a suspended state such that it's only
# waiting for callbacks (including the OnExit one)
# The only way out of this state is to call Exit()
#
OnExit(cbQuit); #Call this when the script exits
WaitForExit();
Download