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  " \
							+ " georef Y   Raster Row   Raster Col  Focal median");
	}

# print to file
fprintf(fHandle, "%s\n\n", "Point   Vector X    Vector Y    georef X  " \
							+ " 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", \
					 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", \
					 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