import logging
import os
import plyfile
import numpy
import time
from osgeo import osr, gdal
from laserfarm.utils import check_dir_exists
from laserfarm.pipeline_remote_data import PipelineRemoteData
logger = logging.getLogger(__name__)
[docs]
class GeotiffWriter(PipelineRemoteData):
""" Write specified bands from point cloud data into separate geotiff files. """
def __init__(self, input_dir=None, bands=None, label=None):
self.pipeline = ('parse_point_cloud',
'data_split',
'create_subregion_geotiffs')
self.InputTiles = []
self.subtilelists = []
self.LengthDataRecord = 0
self.xResolution = 0
self.yResolution = 0
if input_dir is not None:
self.input_path = input_dir
if bands is not None:
self.bands = [bands] if isinstance(bands, str) else bands
if label is not None:
self.label = label
[docs]
def parse_point_cloud(self):
"""
Parse input point cloud and get the following information:
- Tile list
- Length of a single band
- x and y resolution
"""
check_dir_exists(self.input_path, should_exist=True)
# Get list of input tiles
self.InputTiles = [TileFile
for TileFile in os.listdir(self.input_path)
if TileFile.lower().endswith('.ply')]
if not self.InputTiles:
raise IOError('No PLY file in dir: {}'.format(self.input_path))
else:
logger.info('{} PLY files found'.format(len(self.InputTiles)))
# Read one tile and get the template
file = os.path.join(self.input_path, self.InputTiles[0])
template = plyfile.PlyData.read(file)
if not template.elements[0].name == 'vertex':
raise ValueError('Tile PLY file should '
'have vertex as first object')
# Get length of data record (Nr. of elements in each band)
self.LengthDataRecord = len(template.elements[0].data)
logger.info('No. of points per file: {}'.format(self.LengthDataRecord))
# Get resolution, assume a square tile
delta_x = (template.elements[0].data['x'].max()
- template.elements[0].data['x'].min())
delta_y = (template.elements[0].data['y'].max()
- template.elements[0].data['y'].min())
if numpy.isclose(delta_x, 0.) or numpy.isclose(delta_y, 0.):
raise ValueError('Tile should have finite extend in X and Y!')
self.xResolution = (delta_x / (numpy.sqrt(template.elements[0].data['x'].size) - 1))
self.yResolution = (delta_y / (numpy.sqrt(template.elements[0].data['y'].size) - 1))
if not (numpy.isclose(self.xResolution, self.yResolution) and
numpy.isclose(delta_x, delta_y)):
raise ValueError('Tile read is not square!')
logger.info('Resolution: ({}m x {}m)'.format(self.xResolution,
self.yResolution))
return self
[docs]
def data_split(self, xSub, ySub):
"""
Split the input data into sub-regions
:param xSub: number of sub-regions in horizontal direction
:param ySub: number of sub-regions in vertical direction
"""
xcoord = []
ycoord = []
if not self.InputTiles:
raise ValueError('Input tile list is empty!')
for f in self.InputTiles:
comp = f.split('_')
xc = comp[1]
yc = comp[2].split('.')[0]
xcoord.append(xc)
ycoord.append(yc)
# Tile index list
xcint = list(map(float, xcoord))
ycint = list(map(float, ycoord))
# Extent of the tiles
maxxc = max(xcint)
minxc = min(xcint)
maxyc = max(ycint)
minyc = min(ycint)
# Range tile index
xcRange = maxxc - minxc + 1
ycRange = maxyc - minyc + 1
# Range of each sub-region
xcSubRange = numpy.floor(xcRange / xSub)
ycSubRange = numpy.floor(ycRange / ySub)
# Loop per sub-region, find relevant tile of this new tile
# Start from bottom left
logger.info('Splitting data into ({}x{}) sub-regions'.format(xSub,
ySub))
for i in range(xSub):
for j in range(ySub):
if i != xSub - 1 and j != ySub - 1:
# Not the last line/colunm
# Include left/bottom; Exclude right/up
# [Left, Right): [minxc + i*xcSubRange, minxc + (i+1)*xcSubRange);
# [Bottom, top): [minyc + j*ycSubRange, minyc + (j+1)*ycSubRange);
subtiles = [f for k, f in enumerate(self.InputTiles)
if ((minxc + i * xcSubRange) <= xcint[k] < (minxc + (i + 1) * xcSubRange)
and (minyc + j * ycSubRange) <= ycint[k] < (minyc + (j + 1) * ycSubRange))]
if i == xSub - 1 and j == ySub - 1:
# top right corner
# [Left, right]: [minxc + i*xcSubRange; maxxc];
# [Bottom, top]: [minyc + j*ycSubRange; maxyc];
subtiles = [f for k, f in enumerate(self.InputTiles)
if ((minxc + i * xcSubRange) <= xcint[k] <= maxxc
and (minyc + j * ycSubRange) <= ycint[k] <= maxyc)]
if i != xSub - 1 and j == ySub - 1:
# Top line but not top right corner
# [Left, right]: [minxc + i*xcSubRange; minxc + (i+1)*xcSubRange];
# [Bottom, top]: [minyc + j*ycSubRange; maxyc];
subtiles = [f for k, f in enumerate(self.InputTiles)
if ((minxc + i * xcSubRange) <= xcint[k] < (minxc + (i + 1) * xcSubRange)
and (minyc + j * ycSubRange) <= ycint[k] <= maxyc)]
if i == xSub - 1 and j != ySub - 1:
# Right colunm but not top right corner
# [Left, right]: [minxc + i*xcSubRange; maxxc];
# [Bottom, top]: [minyc + j*ycSubRange; minyc + (j+1)*ycSubRange];
subtiles = [f for k, f in enumerate(self.InputTiles)
if ((minxc + i * xcSubRange) <= xcint[k] <= maxxc
and (minyc + j * ycSubRange) <= ycint[k] < (minyc + (j + 1) * ycSubRange))]
self.subtilelists.append(subtiles)
return self
[docs]
def create_subregion_geotiffs(self, output_handle, EPSG=28992):
"""
Export geotiff per sub-region, loop in band dimension
:param output_handle: Handle of output file. The output will be named
as <output_handle>_TILE_<tile ID>_BAND_<band name>
:param EPSG: (Optional) EPSG code of the spatial reference system of
the input data. Default 28992.
"""
check_dir_exists(self.output_folder, should_exist=True)
outfilestem = os.path.join(self.output_folder.as_posix(),
output_handle)
for subTiffNumber in range(len(self.subtilelists)):
infiles = self.subtilelists[subTiffNumber]
logger.info('Processing sub-region GeoTiff no. {} '
'...'.format(subTiffNumber))
logger.info('... number of constituent tiles: '
'{}'.format(len(infiles)))
if infiles:
outfile = '{}_TILE_{:03d}'.format(outfilestem, subTiffNumber)
_make_geotiff_per_band(infiles,
outfile,
self.bands,
self.input_path.as_posix(),
self.LengthDataRecord,
self.xResolution,
self.yResolution,
EPSG)
else:
logger.warning(
'No data in sub-region no. ' + str(subTiffNumber))
logger.info('... processing of sub-region completed.')
return self
def _make_geotiff_per_band(infiles, outfile, band_export, data_directory,
lengthDataRecord, xResolution, yResolution, EPSG):
# Set the coordinate frame
logger.debug('... setting the coordinate frame')
xyData = _plyIntoNumpyArray(data_directory, infiles, lengthDataRecord,
['x', 'y'])
# Shift the coordinates to the center of the cell
xyDataShifted = _shiftTerrain(xyData, xResolution, yResolution)
geoTransform, arrayinfo = _getGeoTransform(xyDataShifted,
xResolution,
yResolution)
# GeoCoding: get the index of each point in the raster
indexX, indexY = _getGeoCoding(xyDataShifted, arrayinfo)
ncols = int(arrayinfo[3])
nrows = int(arrayinfo[7])
for band_name in band_export:
if band_name not in ['x', 'y']:
logger.debug('... creating GeoTiff for band {!s}'.format(band_name))
ct0 = time.time()
# Import one band from PLY
logger.debug('... importing data')
terrainDataOneBand = _plyIntoNumpyArray(data_directory,
infiles,
lengthDataRecord,
[band_name])
# Convert from pointcloud to raster
RasterData = numpy.full((nrows, ncols), numpy.nan)
RasterData[indexY, indexX] = terrainDataOneBand[:, 0]
# Write the single band to geotiff
outfile_band = outfile + "_BAND_" + band_name
_writeGeoTiff(RasterData, band_name, geoTransform, outfile_band,
ncols, nrows, 1, EPSG)
ct1 = time.time()
dct = ct1 - ct0
logger.debug(('... Tiff created in {!s} seconds. Location: '
'{!s}.tif'.format(str(dct), outfile_band)))
def _getGeoTransform(xyData, xres, yres):
"""Adpated to accomodate the orientation expected by geotiffs. """
xmin, ymin, xmax, ymax = [xyData[:, 0].min(), xyData[:, 1].min(),
xyData[:, 0].max(), xyData[:, 1].max()]
ncols = round(((xmax - xmin) / xres) + 1)
nrows = round(((ymax - ymin) / yres) + 1)
geotransform = (xmin, xres, 0, ymax, 0, -1. * yres)
arrayinfo = (xmin, xmax, xres, ncols, ymin, ymax, yres, nrows)
return geotransform, arrayinfo
def _shiftTerrain(terrainData, xres, yres):
"""
This shifts the coordinates by half a cell to account for shift between
target list and cell coordinate assumption made by gdal accommodating
geotiff orientation convention.
"""
tdc = terrainData.copy()
tdx = tdc[:, 0]
tdy = tdc[:, 1]
tdx = tdx - 0.5 * xres
tdy = tdy - 0.5 * yres * (-1.)
tdc[:, 0] = tdx
tdc[:, 1] = tdy
return tdc
def _getGeoCoding(xyData, arrayinfo):
"""Geocoding the point-wise x/y to a raster grid. """
x_idx = (xyData[:, 0] - arrayinfo[0]) / arrayinfo[2]
y_idx = -(xyData[:, 1] - arrayinfo[5]) / arrayinfo[6]
assert numpy.allclose(x_idx, numpy.rint(x_idx)), 'Geo coding failed!'
assert numpy.allclose(y_idx, numpy.rint(y_idx)), 'Geo coding failed!'
return x_idx.astype(int).tolist(), y_idx.astype(int).tolist()
def _writeGeoTiff(featureArrays, bandName, geoTransform, outputFileName, ncols,
nrows, nbands, EPSG_code):
# TODO: READ EPSG_code FROM INPUT PLY
output_raster = gdal.GetDriverByName('GTiff').Create(
outputFileName+".tif", ncols, nrows, nbands, gdal.GDT_Float32,
['COMPRESS=LZW'])
output_raster.SetGeoTransform(geoTransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(EPSG_code)
output_raster.SetProjection(srs.ExportToWkt())
output_raster.SetMetadata({'band': bandName})
rb = output_raster.GetRasterBand(1)
rb.SetMetadata({"band_key": bandName})
rb.WriteArray(featureArrays)
output_raster.FlushCache()
def _plyIntoNumpyArray(directory, tileList, gridLength, columnList):
terrainData = numpy.empty((gridLength * len(tileList), len(columnList)))
for i, file in enumerate(tileList):
if i % 25 == 0 or i == len(tileList) - 1: # first, every 25, and last
logger.debug('... processing tile ' + str(i+1) + ' of ' + str(len(tileList)))
plydata = plyfile.PlyData.read(directory + "/" + file)
for j, column in enumerate(columnList):
terrainData[gridLength * i:gridLength * i + gridLength, j] = plydata.elements[0].data[column]
return terrainData
