LASextractByRegion.sml

# LASextractByRegion.sml

# sample script demonstrates extracting points from an LAS LIDAR file
# using a region and writing to a new LAS file.

# requires TNTmips version 2010 or higher.
# version 16 December 2009

numeric i, j;
class RVC_SHAPE lasIn;					# the input LAS file linked as a shape object
class RVC_OBJITEM objItemIn;
class SR_COORDREFSYS crs;			# Coordinate reference system of input LAS file
class RECT3D lasExtents;				# extents of the input LAS file in its CRS
class REGION Reg;								# region object selected

###########################
### procedure to select the region object and check that its extents overlap those of the LAS file
proc selectRegion () {
	GetInputRegion(Reg);
	Reg.ConvertTo(crs);			# convert to CRS of LAS file

	if (lasExtents.Overlaps(Reg.Extents) == 0) 	# check for overlap; if not, call procedure again to reselect region
		{
		PopupMessage("Region selected does not overlap LAS file extents; please select another ");
		selectRegion();
		}
	}		# end selectRegion()


#####################   MAIN PROGRAM   #############################
clear();

# get input LAS shape object
DlgGetObject("Select input LAS shape object:", "Shape", objItemIn, "ExistingOnly");
lasIn.Open(objItemIn, "Read");

class RVC_GEOREFERENCE georef;		# get default georeference from input LAS file
lasIn.GetDefaultGeoref(georef);

crs = georef.GetCoordRefSys();
printf("Input CRS: %s\n", crs.Name);

# get the extents of the LAS shape object for comparison with region
lasIn.GetExtents(lasExtents);

# open input shape database and main table (table number = 0) for read
class RVC_DBASE_SHAPE dbIn;
dbIn.OpenAsSubobject(lasIn, "Read");

class RVC_DBTABLE tableIn;
tableIn.Open(dbIn, 0, "Read");
printf("Number of Fields in input LAS file = %d\n", tableIn.GetNumFields() );

# call user-defined procedure to select and check the extraction region
selectRegion();

# get filepath for output LAS file
class FILEPATH path = GetOutputFileName("output.las", "Select LAS file to make:", "las");

# make output LAS file for extracted points; use method that takes the RVC_DBTABLE class
# instance for the existing LAS file to set the same point data record type for the new LAS file
class RVC_SHAPE lasOut;
lasOut.MakeLAS(path, crs, tableIn);

class RVC_GEOREFERENCE georefOut;
lasOut.GetDefaultGeoref(georefOut);
printf("Input CRS: %s\n", georefOut.GetCoordRefSys().Name);
 
# open shape database and main table for write
class RVC_DBASE_SHAPE dbOut;
dbOut.OpenAsSubobject(lasOut, "Write");

class RVC_DBTABLE tableOut;
tableOut.Open(dbOut, 0, "Write");
printf("Number of Fields in output LAS file = %d\n", tableOut.GetNumFields() );

# record class instances for reading, copying, and writing records
class RVC_DBTABLE_RECORD recordIn(tableIn);
class RVC_DBTABLE_RECORD recordOut(tableOut);
class RVC_RECORDNUM recordNum;		# container for record number

class STATUSCONTEXT status;
class STATUSDIALOG statusDLG;
statusDLG.Create();
status = statusDLG.CreateContext();
status.BarInit(tableIn.GetNumRecords(), 0);
status.Message = "Processing LIDAR points...";

class POINT2D pt;

# loop through the LIDAR point records to find points inside region
for i = 1 to tableIn.GetNumRecords()
	{
	if (i % 100000 == 0) then
	 	printf("Processing record %d of %d\n", i, tableIn.GetNumRecords() );
	status.BarUpdate(i, tableIn.GetNumRecords(), 0);

	recordNum.Number = i;
	tableIn.Read(recordNum, recordIn);		# read record from input LAS to record container

	pt.x = recordIn.GetValue("X");		# get map position of current point
	pt.y = recordIn.GetValue("Y");

	if (Reg.IsPointInside(pt) )		# check that these map coordinates are inside the extraction region
		{
		recordIn.CopyTo(recordOut);			# copy field values from record in input to a new record for the output

		tableOut.AddRecord(recordOut);		# write the new record to output LAS file
		}
	}

statusDLG.Destroy();

printf("Number of points transferred to output = %d\n", tableOut.GetNumRecords() );
print("Done");