vecfocal.sml

# VecFocal.sml
#
# for each corresponding point in a vector layer - calculates focal mean around that
# same point in a corresponding raster layer
#
# a sample case would be extract maximum points from an elevation raster then use 
# the resulting vector points and the original raster
#
# a sample RVC project file called VecFocal.rvc has these two such files
#
#
# AUTHOR: 			David Wilson
# 						MicroImages, Inc.
# REQUESTED BY:
# CREATION DATE: 	Jan. 20, 1997
# REVISION LIST: 
#
# procedure and function definitions
#
#
# Global Variables
#
raster R;
vector V;
# set to true if you want to see output on console
numeric debugPrintToConsole = true;
# get the input raster and vector
GetInputRaster(R);
GetInputVector(V);
# get the desired output file name
string outFileName$ = PopupString("Input output file Name: ");
# open file for write
class FILE fHandle;
fHandle = fopen(outFileName$, "w");
clear(); # clear the console
# get georef object id's for raster and vector
class Georef georefR, georefV;
georefR = GetLastUsedGeorefObject(R);
georefV = GetLastUsedGeorefObject(V);
# determine number of points in vector
numeric numVectorPoints = NumVectorPoints(V);
array numeric focalOut[numVectorPoints];	# array to output for each point
# print to screen
if (debugPrintToConsole) {
	printf("%s\n\n", "Point   Vector X    Vector Y    georef X  " $2
							+ " georef Y   Raster Row   Raster Col  Focal median");
	}
# print to file
fprintf(fHandle, "%s\n\n", "Point   Vector X    Vector Y    georef X  " $2
							+ " georef Y   Raster Row   Raster Col  Focal median");
# loop over all points - compute mean in 5x5 area around each point
numeric i;
for i = 1 to numVectorPoints {
	# get the vector coordinates
	numeric x, y, xVector, yVector;
	x = V.point[i].Internal.x;
	y = V.point[i].Internal.y;
	# convert vector coordinates to georef
	ObjectToMap(V, x, y, georefV, xVector, yVector);
	# convert raster  coordinates to georef
	# note the reversal of rCol and rLine
	numeric rLine, rCol;
	MapToObject(georefR, xVector, yVector, R, rCol, rLine);
	# round rCol, rLine - could use floor() to force to upper left corner
	rLine = round(rLine);
	rCol = round(rCol);
	# calculate focal mean
	focalOut[i] = FocalMean(R, 5, 5, rLine, rCol);
	# print to screen
	if (debugPrintToConsole) {
		printf("%5d %10.2f  %10.2f  %10.2f  %10.2f  %10d   %10d %e\n", $2
					 i, x, y, xVector, yVector, rLine, rCol, focalOut[i]);
		}
	# print to file
	fprintf(fHandle, "%5d %10.2f  %10.2f  %10.2f  %10.2f  %10d   %10d %e\n", $2
					 i, x, y, xVector, yVector, rLine, rCol, focalOut[i]);
	}
GeorefFree(georefR);	# must free up id's
GeorefFree(georefV);
CloseRaster(R);	# close raster
fclose(fHandle);	# close text file