LAS_to_Raster.sml

# LasToRaster.sml
#
# This script takes one or more LAS files as input and creates Min,
# Max and Count rasters from it.
#
# Method:
# -------
# 
# 1. Calculates total extents of all input Las files and creates 
#    georeference for the output rasters using the same CRS as the
#    Las files.
# 2. Loops through each Las file and updates values in the three 
#    rasters as required.
#
#
# Terranean Mapping Technologies, Brisbane, Australia: 
# 12/5/2009
#####################################################################
clear();
class GEOREF georef;
class RASTER Min,Max,Count;
class RVC_DBTABLE lastable[];
class RVC_SHAPE lasshape[];
class RVC_DBTABLE_RECORD lasRecord;
class RVC_DESCRIPTOR descriptor;
class RVC_OBJITEM lasfiles[];
class STRINGLIST laslist, rastlist;
class RVC_DBASE_ELEMENT dbase;
class RVC_GEOREFERENCE lasref;
class POLYLINE polyline;
class RECT3D extents;
class POINT3D p;
numeric x,y,z;
numeric col,line,l,i;
numeric minX, maxX, minY, maxY, cols, lines, numpts;
array numeric xsrc[4],ysrc[4],zsrc[4],xdst[4],ydst[4],zdst[4];
string str$, name$, path$, rvc$;
# Select LAS files
DlgGetObjects("Las files", "All",lasfiles, "NewOrExisting");
laslist = lasfiles.GetKeys();
numeric numLas = laslist.GetNumItems();
lastable.clear();
lasshape.clear();
# Get output pixel size
numeric cellsize = PopupNum("Cell Size:",1.0,0.1,100,2);
# Get combined LAS extents.
minX = 999999999;
minY = 999999999;
for l = 1 to numLas
{	lasshape[l].Open(lasfiles[l], "Read");
	descriptor = lasshape[i].GetDescriptor();
	name$ = FileNameGetName(descriptor.GetFullName());
	path$ = FileNameGetPath(descriptor.GetFullName());
	rastlist.AddToEnd(name$);
	lasshape[l].GetExtents(extents);
   lasshape[l].GetDefaultGeoref(lasref);
	if (extents.x1 < minX) minX = extents.x1;
	if (extents.x2 > maxX) maxX = extents.x2;
	if (extents.y1 < minY) minY = extents.y1;
	if (extents.y2 > maxY) maxY = extents.y2;
	if (extents.x2 < minX) minX = extents.x2;
	if (extents.x1 > maxX) maxX = extents.x1;
	if (extents.y2 < minY) minY = extents.y2;
	if (extents.y1 > maxY) maxY = extents.y1;
print(name$);
}
printf("\tExtents: \n\t\t Minimum: %f, %f\n\t\t Maximum: %f, %f\n\n", minX, minY , maxX, maxY);
# Calculate Raster dimensions and  setup georef
cols = ceil((maxX - minX) / cellsize);
lines = ceil((maxY - minY) / cellsize);
maxX = minX + (cols*cellsize);
maxY = minY + (lines*cellsize);
print("min/max x",minX,maxX);
print("min/max y",minY,maxY);
print(lines,"lines.");
print(cols,"columns.");
xsrc[1] = 0;		ysrc[1] = 0;			xdst[1] = minX;	ydst[1] = maxY;
xsrc[2] = cols;	ysrc[2] = 0; 			xdst[2] = maxX;	ydst[2] = maxY;
xsrc[3] = cols;	ysrc[3] = lines; 	xdst[3] = maxX;	ydst[3] = minY;
xsrc[4] = 0; 		ysrc[4] = lines; 	xdst[4] = minX;  ydst[4] = minY;
# Create Rasters
GetOutputRaster(Min,lines,cols,"32-bit float");
rvc$ = GetObjectFileName(Min);
CreateRaster(Max,rvc$,"Max","",lines,cols,"32-bit float");
CreateRaster(Count,rvc$,"Count","",lines,cols,"8-bit unsigned");
georef = CreateControlPointGeoref(Min,4,xsrc,ysrc,zsrc,xdst,ydst,zsrc,zsrc,zsrc,zsrc,zsrc,zsrc,zsrc,lasref.GetCoordRefSys());
CopySubobjects(Min,Max,"georef");
CopySubobjects(Min,Count,"georef");
Min = 9999;
Max = -9999;
#Loop through all selected las files
for l = 1 to  numLas
{	lasshape[l].Open(lasfiles[l], "Read");
	print(l,"of",numLas,name$);
	dbase.OpenAsSubobject(lasfiles[l]);
	#open las table under database
	lastable[l].Open(dbase, 0);
	numpts =  lastable[l].GetNumRecords();
	#Loop through each record in the las database
	for i = 1 to numpts
	{	if (! (i % 10000)) then
			{	str$ = sprintf("point %d of %d points", i,numpts);
				SetStatusMessage(str$);
			}
		lastable[l].Read(i, lasRecord);
		x = lasRecord.GetValue("X");
		y = lasRecord.GetValue("Y");
		z = lasRecord.GetValue("Z");
		#Determine which Line and Column in the output rasters the LiDAR point is on
		col = (x - minX) / cellsize;
		line = (maxY - y) / cellsize;
		if (z < Min[line,col]) Min[line,col] = z;
		if (z > Max[line,col]) Max[line,col] = z;
		Count[line,col] = Count[line,col] + 1;
	}
}
# Finish up
SetNull(Min,9999);
SetNull(Max,-9999);
SetNull(Count,0);
CreatePyramid(Max);	CreateHistogram(Max);
CreatePyramid(Min);	CreateHistogram(Min);
CreatePyramid(Count);	CreateHistogram(Count);
#eof