I'm looking for a way to create a height raster (raster of z values - i.e. not just ground poinits). I can accomplish this with gdal by creating a point vector layer from the numpy array then using RasterizeLayer() with options="BURN_VALUE_FROM=Z". However this takes waaaay too long for las files with 10 million plus points.

I have this that works (taken from this gdal cookbook):

def create_z_raster(las, tif, srs):
    Create a height raster

    las = an open laspy File
    tif = 'path/to/tif'
    srs = string format of spatial reference

    # set up the driver in memory
    driver = ogr.GetDriverByName('Memory')

    # create the data source
    data_source = driver.CreateDataSource('in_mem')  # get in_memory driver

    # create the spatial reference
    sr = osr.SpatialReference()

    # create the layer
    layer = data_source.CreateLayer('points', sr, ogr.wkbPoint)

    # process numpy array of las points
    las_points = np.vstack((las.x, las.y, las.z)).transpose()

    print("creating point vectors...")
    start_time = time.time()
    for pnt in las_points:
        # create the feature
        feature = ogr.Feature(layer.GetLayerDefn())
        # create the wkt for the feature
        wkt = 'POINT ({0} {1} {2})'.format(float(pnt[0]), float(pnt[1]), float(pnt[2]))
        # create the point from the wkt
        point = ogr.CreateGeometryFromWkt(wkt)
        # set the feature geometry using the point
        # create the feature in the layer
        # dereference the feature
        feature = None
    print(time.time() - start_time)

    # calculate raster resolutions
    pixel_size = .5
    xmin, xmax, ymin, ymax = layer.GetExtent()
    x_res = int((xmax - xmin) / pixel_size) + 1  # round up and add additional pixel for remainder
    y_res = int((ymax - ymin) / pixel_size) + 1  # round up and add additional pixel for remainder

    # create destination source
    target_ds = gdal.GetDriverByName('GTiff').Create(tif, x_res, y_res, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform((xmin, pixel_size, 0, ymax, 0, -pixel_size))
    band = target_ds.GetRasterBand(1)

    # Rasterize
    gdal.RasterizeLayer(target_ds, [1], layer, options=["BURN_VALUE_FROM=Z"])
    layer = None
    target_ds = None

I would like to know how to do this without first creating point vectors. How do I build a numpy array of z values that can be written to a raster band with a width of x_res and a height of y_res from a las file with n number of points? Or is this a wildly inappropriate approach?

  • 1
    So are you just looking to interpolate a LAS file to a raster? – WhiteboxDev May 31 '18 at 20:35
  • It's a valid approach except to minimize overhead you should be resolving the record to a cell rather than generating a point; start with a numpy array the size of your raster, with origin and cell size, read each LAS record and calculate the cell it falls on, test that index in your array and populate or overwrite then when you're all done reading create a raster and write your array to it. I have written code to do just what you are trying to achieve using the GDAL libs but coding in C++ is the only way to get any sort of reasonable processing time; can you code C++? – Michael Stimson Jun 1 '18 at 4:10
  • Are you interested in using a CLI? I use FUSION or the lidr package in R to calculate these metrics. – Aaron Jul 6 '18 at 13:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.