Syntax Highlighing:
comments, key words, predefined symbols, class members & methods, functions & classes
# SOILTEST.SML
#
# SML script to perform surface fitting on a set of attribute values
# attached to vector points that represent an array of soil sample locations.
# Values are read from a attached database table called "TABLE".
# The script uses a field boundary polygon in a second vector object to make a mask raster.
# The polygon is specified by element number, with the default being 1.
# The mask raster is applied to each output surface raster to assign null for cells outside
# the field boundary. The mask raster is also used as a reference raster for the surface
# fitting function to ensure that the surface rasters cover the full map extents
# of the field. Written for TNTmips v6.30.
# Author: Randy Smith, MicroImages, Inc.
# Date: 10 March 2000
###########################################################################################
# Revisions
###########################################################################################
# 15 March 2000: MakeSurf procedure modified to delete raster pyramids automatically
# created by surface fitting function prior to applying mask. Add
# popup dialog to prompt for element number of boundary polygon for
# conversion to region to use to make mask.
# 29 July 2003: modified to conform to strict syntax checking. MakeSurf procedure
# modified to take valuePoint$ parameter instead of using global
# parameter value.
#
# 2 December 2003: modified to set 0 as null for mask (mask for some reason is no longer
# automatically created with 0 set as null).
#
# 21 September 2004: modified for v7.0 to use internal REGION2D instead of region object
# to set up the mask raster.
###########################################################################################
# global variable declarations
###########################################################################################
vector Boundary, Points;
raster Surface, Mask;
class Georef vecGeoref;
string surfName$;
string surfDescript$;
string rasttype$;
string ptQuery$;
string filename$;
string filedescript$;
numeric cellsize, horizCell, vertCell;
numeric boundElemNum;
###########################################################################################
# Assign fixed values for global variables for Minimum Curvature surface fitting
###########################################################################################
numeric usePoints = 1;
numeric useLines = 0;
string searchMethod$ = "Square";
numeric searchDist = 50;
numeric weightPower = 2.0;
numeric tension = 0.5;
string selectPoints$ = "1";
string selectLines$ = "";
string valueLine$ = "";
###########################################################################################
# Define procedures
###########################################################################################
######################################
# Procedure to make surface raster
######################################
proc MakeSurf(string name$, string descript$, string type$, string valuePoints$) {
# Create blank raster for surface fitting with desired file / object names,
# data type, and cell size.
CreateRasterFromObject(Boundary, Surface, filename$, name$, descript$,
horizCell, vertCell, type$, vecGeoref);
Surface = 0;
CloseVector(Boundary);
# perform surface fitting from vector Points to raster Surface
# using mask raster as a reference raster to control extents and
# georeferencing
SurfaceFitMinimumCurvature(Points, Surface, usePoints, useLines,
searchMethod$, searchDist, weightPower, tension, selectPoints$,
selectLines$, valuePoints$, valueLine$, Mask);
DeletePyramid(Surface); # Delete pyramids created by surface fitting function
# apply mask to raster Surface
for each Surface {
Surface = Surface * Mask;
}
SetNull(Surface, 0);
CreateHistogram(Surface);
CreatePyramid(Surface, 1);
CloseRaster(Surface);
}
#############################################################################################
# Get input objects
#############################################################################################
### Prompt for input object selection
PopupMessage("Select a vector object containing soil test points \n and one containing the field boundary polygon");
GetInputVector(Points); # vector point object
GetInputVector(Boundary); # vector object containing field boundary polygon
vecGeoref = GetLastUsedGeorefObject(Boundary);
boundElemNum = PopupNum("Enter element number of boundary polygon", 1, 1, 10000, 0);
# Create Project File to contain output surface rasters
filename$ = GetOutputFileName("","Enter name of output project file:", "rvc");
filedescript$ = PopupString("Enter description of output project file","Surface rasters from attributes of soil test points");
CreateProjectFile(filename$,filedescript$);
# Prompt for desired cell size of output surface rasters
cellsize = PopupNum("Desired cell size of surfaces in meters:",2,0,100,0);
horizCell = cellsize;
vertCell = cellsize;
#############################################################################################
# Make temporary in-memory region to create mask for surface rasters
#############################################################################################
# Declare region name
class REGION2D RegField; # region from field boundary polygon
# Create region from field boundary polygon
RegField = ConvertVectorPolyToRegion(Boundary,boundElemNum,vecGeoref);
#############################################################################################
# Create a mask raster matching extents of field boundary
# polygon with value = 1 inside field boundary.
#############################################################################################
CreateRasterFromObject(Boundary,Mask,filename$,"Mask","Mask raster with 1 inside field boundary",
horizCell,vertCell,"binary",vecGeoref); # blank raster with desired extents
CloseVector(Boundary);
Mask = 0; # assign initial value to mask
# IgnoreNull(Mask);
for each Mask in RegField {
Mask = 1;
}
#############################################################################################
# Fit surface to point values in TABLE.PH
#############################################################################################
ptQuery$ = "TABLE.PH";
surfName$ = "PH";
surfDescript$ = "PH surface raster";
rasttype$ = "32-bit float";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.BUFFER
#############################################################################################
ptQuery$ = "TABLE.BUFFER";
surfName$ = "BUFFER";
surfDescript$ = "BUFFER surface raster";
rasttype$ = "32-bit float";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.CA_INDEX
#############################################################################################
ptQuery$ = "TABLE.CA_INDEX";
surfName$ = "CA_INDEX";
surfDescript$ = "CA_INDEX surface raster";
rasttype$ = "8-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.PERCA
#############################################################################################
ptQuery$ = "TABLE.PERCA";
surfName$ = "PERCA";
surfDescript$ = "PERCA surface raster";
rasttype$ = "32-bit float";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.CA
#############################################################################################
ptQuery$ = "TABLE.CA";
surfName$ = "CA";
surfDescript$ = "CA surface raster";
rasttype$ = "16-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.MG_INDEX
#############################################################################################
ptQuery$ = "TABLE.MG_INDEX";
surfName$ = "MG_INDEX";
surfDescript$ = "MG_INDEX surface raster";
rasttype$ = "8-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.PERMG
#############################################################################################
ptQuery$ = "TABLE.PERMG";
surfName$ = "PERMG";
surfDescript$ = "PERMG surface raster";
rasttype$ = "32-bit float";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.MG
#############################################################################################
ptQuery$ = "TABLE.MG";
surfName$ = "MG";
surfDescript$ = "MG surface raster";
rasttype$ = "16-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.K_INDEX
#############################################################################################
ptQuery$ = "TABLE.K_INDEX";
surfName$ = "K_INDEX";
surfDescript$ = "K_INDEX surface raster";
rasttype$ = "8-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.PERK
#############################################################################################
ptQuery$ = "TABLE.PERK";
surfName$ = "PERK";
surfDescript$ = "PERK surface raster";
rasttype$ = "32-bit float";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.K
#############################################################################################
ptQuery$ = "TABLE.K";
surfName$ = "K";
surfDescript$ = "K surface raster";
rasttype$ = "16-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.KMGRATIO
#############################################################################################
ptQuery$ = "TABLE.KMGRATIO";
surfName$ = "KMGRATIO";
surfDescript$ = "K / MG RATIO surface raster";
rasttype$ = "32-bit float";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.K_MG_IND
#############################################################################################
ptQuery$ = "TABLE.K_MG_IND";
surfName$ = "K_MG_IND";
surfDescript$ = "K - MG INDEX surface raster";
rasttype$ = "8-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.P_INDEX
#############################################################################################
ptQuery$ = "TABLE.P_INDEX";
surfName$ = "P_INDEX";
surfDescript$ = "P_INDEX surface raster";
rasttype$ = "8-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.P1
#############################################################################################
ptQuery$ = "TABLE.P1";
surfName$ = "P1";
surfDescript$ = "P1 surface raster";
rasttype$ = "8-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.BICARB
#############################################################################################
ptQuery$ = "TABLE.BICARB";
surfName$ = "BICARB";
surfDescript$ = "BICARB surface raster";
rasttype$ = "8-bit unsigned";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.CEC
#############################################################################################
ptQuery$ = "TABLE.CEC";
surfName$ = "CEC";
surfDescript$ = "CEC surface raster";
rasttype$ = "32-bit float";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
#############################################################################################
# Fit surface to point values in TABLE.OM
#############################################################################################
ptQuery$ = "TABLE.OM";
surfName$ = "OM";
surfDescript$ = "Organic matter surface raster";
rasttype$ = "32-bit float";
MakeSurf(surfName$, surfDescript$, rasttype$, ptQuery$);
# End
Download