VIEWSHED.sml

# VIEWSHED.sml
# Computes binary viewshed raster from elevation raster input
# and vector object containing lines.

# Computes viewshed for all vertices on all lines
# of a corresponding Vector object
# Could use to get viewshed for all points
# on one or more roads.

### Declarations
array numeric xArray[10];	 # use to hold line vertices
array numeric yArray[10];

class RASTER Rin, Rout;
class VECTOR V;

class GEOREF georefR, georefV;
class TRANS2D_MAPGEN vecGeoToRastGeo; 
numeric numLines, numPoints;
numeric thisLine, numPointsInThisLine;
class MATRIX matX, matY;
numeric currentPoint, i;
numeric xVector, yVector, rCol, rLine;
class POINT2D vecMapCoord, rastMapCoord;
numeric percent, height, zScale;

####
clear(); 	# clear the console

# get the input and output raster and Vector object
GetInputRaster( Rin );
GetInputVector( V );
GetOutputRaster( Rout, NumLins( Rin ), NumCols( Rin ), "binary" );

# get georef object id's for raster and vector
georefR = GetLastUsedGeorefObject( Rin );
georefV = GetLastUsedGeorefObject( V );
vecGeoToRastGeo.InputProjection = georefV.Projection;
vecGeoToRastGeo.OutputProjection = georefR.Projection;

# first count total number of vertices in all lines
numLines = NumVectorLines( V );
numPoints = 0;

for thisLine = 1 to numLines {
	numPointsInThisLine = GetVectorLinePointList( V, xArray, yArray, thisLine );
	numPoints += numPointsInThisLine;
	} 

# now create matrices large enough to hold points
matX = CreateMatrix(1, numPoints );
matY = CreateMatrix(1, numPoints );

# now loop through lines and fill in matrices with points

currentPoint = 0; 	# this will be our matrix index

for thisLine = 1 to numLines {
	numPointsInThisLine = GetVectorLinePointList( V, xArray, yArray, thisLine );

	for i = 1 to numPointsInThisLine {
		# convert vector object coordinates to map coordinates from its georef
		ObjectToMap( V, xArray[i], yArray[i], georefV, xVector, yVector );

		# convert vector map coordinates to map coordinates in raster's georef
		vecMapCoord.x = xVector;		# first assign map coordinates to a Point class
		vecMapCoord.y = yVector;		# required by TRANS2D_MAPGEN

		# do conversion using method on class TRANS2D_MAPGEN
		rastMapCoord = vecGeoToRastGeo.ConvertPoint2DFwd(vecMapCoord);

		# convert vector coordinates to raster coordinates
		# note the order of rCol and rLine
		# col = x coord, row = y coord
		MapToObject( georefR, rastMapCoord.x, rastMapCoord.y, Rin, rCol, rLine );

		# force to upper left corner
		rLine = floor( rLine ); 
		rCol = floor( rCol );

		# want matrix to index from zero to numPoints
		SetMatrixItem( matX, 0, currentPoint, rCol );
		SetMatrixItem( matY, 0, currentPoint, rLine );
		currentPoint += 1;
		} # end of for i

	} # end of for each line in V

# set default values - want any one point to make "visible"
percent = 100 / numPoints;
height = 2;
zScale = 1;

# now create the viewshed Raster object
RasterToBinaryViewshed(  Rin, Rout, numPoints, matX, matY, percent, height, zScale );


$ifdef _MI_INTERNAL_AUTOTEST_MODE  # This is used for auto testing purposes only
	numeric checksum = 0;
	for each Rout {
		if (!IsNull(Rout)) checksum = checksum + Rout;
		}
	print("checksum:",checksum);
$endif

# clean up
GeorefFree( georefR );	 # must free up id's
GeorefFree( georefV );
DestroyMatrix( matX );
DestroyMatrix( matY );
CloseRaster( Rin );
CloseRaster( Rout );