SOILTEST.sml

#	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