No values in output when applying processor directly to a LAS file.

[gillins]

Original report by Daniel Clewley (Bitbucket: danclewley, GitHub: danclewley).

---

I think there is a problem applying the LiDAR processor directly to LAS files as the processor example from http://pylidar.org/en/latest/processorexamples.html produces a blank image when applied directly to a LAS file but the expected result when the file is converted to SPD format first.

I've tested on a couple of LAS files I have and apl1dr_x509000ys6945000z56_2009_ba1m6_pbrisba.las from the pylidar test data and get the same result. For all files I have generated a spatial index first using lasindex. I'm using pylidar 0.3.0 installed via conda on Python 3.5 under Linux.

It looks like the problem is getPointsByBins always returns 0 points for the LAS file but not the SPD version. The code below will reproduce the problem by running the same function on a LAS file and an SPD v4 created from the input LAS. Adding a line to print zValues.shape shows no points being found for the LAS file.

import os
import numpy
from rios import cuiprogress
from pylidar import lidarprocessor
from pylidar.toolbox.translate.las2spdv4 import translate
from pylidar.lidarformats import generic

def writeImageFunc(data):
     zValues = data.input_points.getPointsByBins(colNames='Z')
     (maxPts, nRows, nCols) = zValues.shape
     nullval = 0
     if maxPts > 0:
         minZ = zValues.min(axis=0)
         stack = numpy.ma.expand_dims(minZ, axis=0)
     else:
         stack = numpy.empty((1, nRows, nCols), dtype=zValues.dtype)
         stack.fill(nullval)
     data.output_image.setData(stack)

if __name__ == '__main__':

     input_las = 'testdata/apl1dr_x509000ys6945000z56_2009_ba1m6_pbrisba.las'
     imported_spd = os.path.splitext(input_las)[0] + '.spd'

     info = generic.getLidarFileInfo(input_las)

    importedSPD = imported_spd
    translate(info, input_las, imported_spd, epsg=27700,
              pulseIndex='FIRST_RETURN', spatial=True, binSize=2.0, 
              buildPulses=True)

    for in_points in [input_las, imported_spd]:
        out_img = in_points + '.img'

        dataFiles = lidarprocessor.DataFiles()
        dataFiles.input_points = lidarprocessor.LidarFile(in_points,
                                                         lidarprocessor.READ)
        dataFiles.output_image = lidarprocessor.ImageFile(out_img,
                                                          lidarprocessor.CREATE)

        if os.path.splitext(in_points)[-1].lower() != '.spd':
            dataFiles.input_points.setLiDARDriverOption('BIN_SIZE', 2.0)

        controls = lidarprocessor.Controls()
        controls.setProgress(cuiprogress.CUIProgressBar())

        lidarprocessor.doProcessing(writeImageFunc, dataFiles)

[gillins]

Original comment by Sam Gillingham (Bitbucket: gillins, GitHub: gillins).

---

Hi Dan,

Thanks for the report. We are thinking of removing spatial processing support (see https://groups.google.com/forum/?utm_medium=email&utm_source=footer#!msg/pylidar/9n4sp_G2r3c/DiceptZdEwAJ) so I'm probably not going to be able to justing spending any time on fixing this. Feel free to see if you can fix it in the interim if it is really bothering you, otherwise use non-spatial mode and hold the whole image in memory updating cells as needed (the new approach we would recommend).

Please contribute to the google group discussion if you have any thoughts on spatial processing.

Thanks, Sam.

[gillins]

Original comment by Daniel Clewley (Bitbucket: danclewley, GitHub: danclewley).

---

Thanks Sam.

I've tried setting controls.setSpatialProcessing(False) and I get the error Can only process image inputs when doing spatial processing. Is it possible to do non-spatial processing and have an image output? I couldn't figure out how do this.

I saw #11 had been opened about spatial processing. Once I've had more of chance to understand pylidar I'll add some of my thoughts on this.

[gillins]

Original comment by Sam Gillingham (Bitbucket: gillins, GitHub: gillins).

---

Hi Dan,

Under the proposed new system you won't be able to read/write images using pylidar. I have adapted @neilflood 's example script below to show the steps you would have to go through to save the image yourself. It is likely we would write some helper functions so not so much code needs to be written each time.

Any thoughts appreciated.

Sam.

#!python

def getCmdargs():
    """
    Get commandline arguments
    """
    p = argparse.ArgumentParser()
    p.add_argument("infile", help=("Input SPDV4 file."))
    p.add_argument("-r", "--resolution", default=0.1, type=float,
        help="Resolution (i.e. pixel size) (in metres) of DEM image to work with (default=%(default)s)")
    p.add_argument("--blocksize", default=256, type=int,
        help=("Pylidar processing blocksize (default=%(default)s). For bizarre reasons, pylidar "+
            "uses the square of this as the number of pulses to work with on each block"))
    p.add_argument("--xmin", type=float,
        help="X coord of western edge of DEM image. Default is X_MIN of input file. ")
    p.add_argument("--xmax", type=float,
        help="X coord of eastern edge of DEM image. Default is X_MAX of input file. ")
    p.add_argument("--ymin", type=float,
        help="Y coord of southern edge of DEM image. Default is Y_MIN of input file. ")
    p.add_argument("--ymax", type=float,
        help="Y coord of northern edge of DEM image. Default is Y_MAX of input file. ")
    p.add_argument("--saveminz", 
        help="Img file in which to save minz image (default will not save this)")

def createInitialMinzImg(cmdargs, zNull):
    """
    Create an image array of the minimum Z value for each pixel. Return a 
    2-d array of minZ values for each pixel. 
    """
    dataFiles = lidarprocessor.DataFiles()

    dataFiles.input1 = lidarprocessor.LidarFile(cmdargs.infile, lidarprocessor.READ)  

    otherargs = lidarprocessor.OtherArgs()

    controls = lidarprocessor.Controls()
    controls.setWindowSize(cmdargs.blocksize)
    controls.setSpatialProcessing(False)

    # Set up the minZ image array, along with appropriate coordinates relating it to 
    # the real world (i.e. xMin, yMax, res). 
    otherargs.res = cmdargs.resolution
    otherargs.zNull = zNull
    otherargs.xMin = cmdargs.xmin
    otherargs.yMax = cmdargs.ymax
    nCols = int(numpy.ceil((cmdargs.xmax - cmdargs.xmin) / otherargs.res))
    nRows = int(numpy.ceil((cmdargs.ymax - cmdargs.ymin) / otherargs.res))
    otherargs.minZ = numpy.zeros((nRows, nCols), dtype=numpy.float32)
    otherargs.minZ.fill(otherargs.zNull)

    lidarprocessor.doProcessing(doMinZ, dataFiles, otherArgs=otherargs, controls=controls)

    return otherargs.minZ

def doMinZ(data, otherargs):
    x = data.input1.getPoints(colNames='X')
    y = data.input1.getPoints(colNames='Y')
    z = data.input1.getPoints(colNames='Z')

    updateMinZ(x, y, z, otherargs.minZ, otherargs.xMin, 
        otherargs.yMax, otherargs.res, otherargs.zNull)

@jit
def updateMinZ(x, y, z, minzImg, xMin, yMax, res, zNull):
    """
    Called from doProcessing(). 

    For the current block of points, given by the x, y, z coord arrays, update the minZ
    image array. Note that the minZ array starts full of nulls, and each time this function 
    is called (for the next block of points), we update whichever pixels are relevant. 
    """
    numPts = x.shape[0]
    (nRows, nCols) = minzImg.shape

    # Loop explicitly over all points. 
    for i in range(numPts):
        # Calculate the row/column in the image, for the (x, y) coords of the current point
        r = int((yMax - y[i]) / res)
        c = int((x[i] - xMin) / res)

        if r >= 0 and r < nRows and c >= 0 and c < nCols:
            # If this is the first point in this pixel, or if this new point has lower Z than
            # previous, then replace the minZ value
            if minzImg[r, c] == zNull or minzImg[r, c] > z[i]:
                minzImg[r, c] = z[i]

def main():
    """
    Main routine
    """
    cmdargs = getCmdargs()
    info = generic.getLidarFileInfo(cmdargs.infile).header
    defaultBounds(cmdargs, info)

    zNull = -10000

    print("Creating minZ")
    minZ = createInitialMinzImg(cmdargs, zNull)
    if cmdargs.saveminz is not None:
        saveImage(cmdargs.saveminz, minZ, cmdargs.xmin, 
            cmdargs.ymax, cmdargs.resolution, info['SPATIAL_REFERENCE'], zNull)

def saveImage(outfile, imgArr, xmin, ymax, res, proj, nullVal):
    """
    Save the given image array as a GDAL raster file. Only saves a single band. 
    """
    (nRows, nCols) = imgArr.shape

    gdalTypeDict = {numpy.dtype(numpy.float32):gdal.GDT_Float32, 
        numpy.dtype(numpy.uint8):gdal.GDT_Byte}
    gdalType = gdalTypeDict[imgArr.dtype]

    drvr = gdal.GetDriverByName('HFA')
    ds = drvr.Create(outfile, nCols, nRows, 1, gdalType, ['COMPRESS=YES'])
    band = ds.GetRasterBand(1)
    band.WriteArray(imgArr)
    band.SetNoDataValue(nullVal)

    ds.SetProjection(proj)
    ds.SetGeoTransform((xmin, res, 0, ymax, 0, -res))

def defaultBounds(cmdargs, info):
    """
    Fill in default bounds for image
    """
    # Default values for xmin and ymax
    if cmdargs.xmin is None or cmdargs.xmax is None or cmdargs.ymin is None or cmdargs.ymax is None:
        if cmdargs.xmin is None:
            cmdargs.xmin = info['X_MIN']
        if cmdargs.xmax is None:
            cmdargs.xmax = info['X_MAX']
        if cmdargs.ymin is None:
            cmdargs.ymin = info['Y_MIN']
        if cmdargs.ymax is None:
            cmdargs.ymax = info['Y_MAX']

if __name__ == "__main__":
    main()