LiDAR Peak Analysis: What It Takes

[Boggy B]

Waking up this thread again to post something I recently found useful in lidar analysis. In the past I've used QGIS for analysis and verification, but lately I was hunting for something to explore lidar in 3D without all the fiddlyness and performance issues of QGIS. I stumbled on displaz (https://c42f.github.io/displaz/), a super barebones lidar viewer that nonetheless provides enough functionality out of the box for rapid, more spatially-oriented analysis.

I was looking at some unranked tower on Sneffels and was shocked at how easy it is to pinpoint the summit/saddle points with this tool. Load up the LAZ (from TNM), middle-click to select and center rotation on any point, left-click+drag to rotate, right-click+drag or scrollwheel to zoom (Windows). Here's a view of the tower I was looking at ("trim radius" under Shader Parameters limits rendering around the selected point but it's mostly cosmetic):

displaz1.jpg

Rotating/zooming around any point near the top made it obvious there were only a couple summit candidates. In the Log to right you can see attributes of the selected point. The format of "position" is X (lon UTM meters) Y (lat UTM meters) Z (height in meters). "vector diff" shows the delta from the previously-selected point in the same format.

displaz2.jpg

Saddles are also easy to home in. On this particular feature, there's a very sharp ridge in the saddle, and it's clear from rotating around that lidar didn't hit the exact lowest point. I'm not sure how you all are deciding which one to use in this case.

displaz3.jpg

The tool converts the X Y values from the source CRS (NAD83/Conus Albers in this case) to UTM, so to get recognizable WGS84 coordinates you'll have to convert them. There are options detailed earlier in this thread but an easy one is to enter them to "Standard UTM" here (http://rcn.montana.edu/resources/Converter.aspx), 13N will be the zone for most of CO except the western edge (west of Durango).

This method is probably horrendous if the summit/saddle is a large flattish area, though you can edit the shader code to colorize by the Z value. You can also natively colorize the points by classification, with all the caveats mentioned upthread. I wonder if the 3D visualization would help with identifying high points under canopy though.

[bdloftin77]

Waking up this thread again to post something I recently found useful in lidar analysis. In the past I've used QGIS for analysis and verification, but lately I was hunting for something to explore lidar in 3D without all the fiddlyness and performance issues of QGIS. I stumbled on displaz (https://c42f.github.io/displaz/), a super barebones lidar viewer that nonetheless provides enough functionality out of the box for rapid, more spatially-oriented analysis.

Thanks for sharing!

[Boggy B]

Below is shader code (just modified the default one) which gives you the ability to define height-based gradients, and, more usefully, upper/lower bounds adjustable to the cm in realtime.
Save as .glsl and Shader > Open, or dump into Shader > Edit and then apply with Shader > Compile.
Here's what it looks like analyzing a peak with the bounds adjusted to isolate the summit/saddle points (white at top is the highpoint):

shader.jpg

Code: Select all

#version 150

uniform mat4 modelViewMatrix;
uniform mat4 projectionMatrix;
uniform mat4 modelViewProjectionMatrix;

//------------------------------------------------------------------------------
#if defined(VERTEX_SHADER)

uniform float pointRadius = 0.1;    //# uiname=Point Radius; min=0.001; max=10
uniform float trimRadius = 1000000; //# uiname=Trim Radius; min=1; max=1000000
uniform float reference = 400.0;    //# uiname=Reference Intensity; min=0.001; max=100000
uniform float exposure = 1.0;       //# uiname=Exposure; min=0.001; max=10000
uniform float contrast = 1.0;       //# uiname=Contrast; min=0.001; max=10000
uniform int colorMode = 0;          //# uiname=Color Mode; enum=Intensity|Colour|Return Index|Point Source|Las Classification|File Number|Z Value
uniform int minZ = 0;               //# uiname=Min Z (m); min=0; max=1000000
uniform int minZFract = 0;          //# uiname=Min Z +/- (cm); min=-99; max=99
uniform int maxZ = 1000000;         //# uiname=Max Z (m); min=0; max=1000000
uniform int maxZFract = 0;          //# uiname=Max Z +/- (cm); min=-99; max=99
uniform int gradientMin = 6;        //# uiname=Gradient Min; enum=White|Black|Red|Orange|Yellow|Green|Blue|Purple
uniform int gradientMid = 0;        //# uiname=Gradient Mid; enum=None|White|Black|Red|Orange|Yellow|Green|Blue|Purple
uniform int gradientMax = 3;        //# uiname=Gradient Max; enum=White|Black|Red|Orange|Yellow|Green|Blue|Purple
uniform int colorBounding = 1;      //# uiname=Color Bounding; enum=None|White|Black|Red|Orange|Yellow|Green|Blue|Purple
uniform int selectionMode = 0;      //# uiname=Selection; enum=All|Classified|First Return|Last Return|First Of Several
uniform float minPointSize = 0;
uniform float maxPointSize = 400.0;

vec3 colors[8];

// Point size multiplier to get from a width in projected coordinates to the
// number of pixels across as required for gl_PointSize
uniform float pointPixelScale = 0;
uniform vec3 cursorPos = vec3(0);
uniform int fileNumber = 0;
in float intensity;
in vec3 position;
in vec3 color;
in int returnNumber;
in int numberOfReturns;
in int pointSourceId;
in int classification;
//in float heightAboveGround;

flat out float modifiedPointRadius;
flat out float pointScreenSize;
flat out float zColorScale;
flat out vec3 pointColor;
flat out int markerShape;
flat out vec3 gradientColor1;
flat out vec3 gradientColor2;
flat out vec3 gradientColor3;
flat out vec3 oobColor;

float tonemap(float x, float reference, float contrast)
{
    float Y = pow(x/reference, contrast);
    Y = Y / (1.0 + Y);
    return Y;
}

vec3 jet_colormap(float x)
{
    if (x < 0.125)
        return vec3(0, 0, 0.5 + 4*x);
    if (x < 0.375)
        return vec3(0, 4*(x-0.125), 1);
    if (x < 0.625)
        return vec3(4*(x-0.375), 1, 1 - 4*(x-0.375));
    if (x < 0.875)
        return vec3(1, 1 - 4*(x-0.625), 0);
    return vec3(1 - 4*(x-0.875), 0, 0);
}

void main()
{
    // init colors
    colors[0] = vec3(1.0, 1.0, 1.0); // White
    colors[1] = vec3(0.0, 0.0, 0.0); // Black
    colors[2] = vec3(1.0, 0.0, 0.0); // Red
    colors[3] = vec3(1.0, 0.5, 0.0); // Orange
    colors[4] = vec3(1.0, 1.0, 0.0); // Yellow
    colors[5] = vec3(0.0, 1.0, 0.0); // Green
    colors[6] = vec3(0.0, 0.0, 1.0); // Blue
    colors[7] = vec3(0.5, 0.0, 0.5); // Purple

    vec4 p = modelViewProjectionMatrix * vec4(position,1.0);
    float r = length(position.xy - cursorPos.xy);
    float trimFalloffLen = min(5, trimRadius/2);
    float trimScale = min(1, (trimRadius - r)/trimFalloffLen);
    modifiedPointRadius = pointRadius * trimScale;
    pointScreenSize = clamp(2*pointPixelScale*modifiedPointRadius / p.w, minPointSize, maxPointSize);
    markerShape = 1;
    // Compute vertex color
    if (colorMode == 0)
        pointColor = tonemap(intensity, reference, contrast) * vec3(1);
    else if (colorMode == 1)
        pointColor = contrast*(exposure*color - vec3(0.5)) + vec3(0.5);
    else if (colorMode == 2)
        pointColor = vec3(0.2*returnNumber*exposure, 0.2*numberOfReturns*exposure, 0);
    else if (colorMode == 3)
    {
        markerShape = (pointSourceId+1) % 5;
        vec3 cols[] = vec3[](vec3(1,1,1), vec3(1,0,0), vec3(0,1,0), vec3(0,0,1),
                             vec3(1,1,0), vec3(1,0,1), vec3(0,1,1));
        pointColor = cols[(pointSourceId+3) % 7];
    }
    else if (colorMode == 4)
    {
        // Colour according to some common classifications defined in the LAS spec
        pointColor = vec3(exposure*classification);
        if (classification == 2)      pointColor = vec3(0.33, 0.18, 0.0); // ground
        else if (classification == 3) pointColor = vec3(0.25, 0.49, 0.0); // low vegetation
        else if (classification == 4) pointColor = vec3(0.36, 0.7,  0.0); // medium vegetation
        else if (classification == 5) pointColor = vec3(0.52, 1.0,  0.0); // high vegetation
        else if (classification == 6) pointColor = vec3(0.8,  0.0,  0.0); // building
        else if (classification == 9) pointColor = vec3(0.0,  0.0,  0.8); // water
    }
    else if (colorMode == 5)
    {
        // Set point colour and marker shape cyclically based on file number
        markerShape = fileNumber % 5;
        pointColor = vec3((1.0/2.0) * (0.5 + (fileNumber % 2)),
                          (1.0/3.0) * (0.5 + (fileNumber % 3)),
                          (1.0/5.0) * (0.5 + (fileNumber % 5)));
    }
    else if (colorMode == 6) // Color mode based on z value
    {		
        // get fractional min/max Z
        float minZReal = minZ + (minZFract / 100.0);
        float maxZReal = maxZ + (maxZFract / 100.0);
	
        float normalizedZ = (position.z - minZReal) / (maxZReal - minZReal); // Normalize z value
		
        // get colors
        if (colorBounding > 0) oobColor = colors[colorBounding - 1];
        gradientColor1 = colors[gradientMin];
        if (gradientMid > 0)
        {
            gradientColor2 = colors[gradientMid - 1];
            gradientColor3 = colors[gradientMax];
            zColorScale = 2.0;
        }
        else
        {
            gradientColor2 = colors[gradientMax];
            zColorScale = 1.0;
        }
		
        if (colorBounding == 0 || (position.z >= minZReal && position.z <= maxZReal))
        {
            if (gradientMid == 0 || normalizedZ < 0.5) pointColor = mix(gradientColor1, gradientColor2, normalizedZ * zColorScale);
            else pointColor = mix(gradientColor2, gradientColor3, (normalizedZ - 0.5) * zColorScale);
            //pointColor = vec3(normalizedZ, 0.5, 1.0 - normalizedZ); // Example color mapping
        }
        else pointColor = oobColor;
    }
    /*
    else if (colorMode == 8)
    {
        // Color based on height above ground
        pointColor = 0.8*jet_colormap(tonemap(0.16*heightAboveGround, exposure, 3.8*contrast));
    }
    */
    if (selectionMode != 0)
    {
        if (selectionMode == 1)
        {
            if (classification == 0)
                markerShape = -1;
        }
        else if (selectionMode == 2)
        {
            if (returnNumber != 1)
                markerShape = -1;
        }
        else if (selectionMode == 3)
        {
            if (returnNumber != numberOfReturns)
                markerShape = -1;
        }
        else if (selectionMode == 4)
        {
            if (returnNumber != 1 || numberOfReturns < 2)
                markerShape = -1;
        }
    }
    // Ensure zero size points are discarded.  The actual minimum point size is
    // hardware and driver dependent, so set the markerShape to discarded for
    // good measure.
    if (pointScreenSize <= 0)
    {
        pointScreenSize = 0;
        markerShape = -1;
    }
    else if (pointScreenSize < 1)
    {
        // Clamp to minimum size of 1 to avoid aliasing with some drivers
        pointScreenSize = 1;
    }
    gl_PointSize = pointScreenSize;
    gl_Position = p;
}


//------------------------------------------------------------------------------
#elif defined(FRAGMENT_SHADER)

uniform float markerWidth = 0.3;

flat in float modifiedPointRadius;
flat in float pointScreenSize;
flat in vec3 pointColor;
flat in int markerShape;

out vec4 fragColor;

// Limit at which the point is rendered as a small square for antialiasing
// rather than using a specific marker shape
const float pointScreenSizeLimit = 2;
const float sqrt2 = 1.414213562;

void main()
{
    if (markerShape < 0) // markerShape == -1: discarded.
        discard;
    // (markerShape == 0: Square shape)
#   ifndef BROKEN_GL_FRAG_COORD
    gl_FragDepth = gl_FragCoord.z;
#   endif
    if (markerShape > 0 && pointScreenSize > pointScreenSizeLimit)
    {
        float w = markerWidth;
        if (pointScreenSize < 2*pointScreenSizeLimit)
        {
            // smoothly turn on the markers as we get close enough to see them
            w = mix(1, w, pointScreenSize/pointScreenSizeLimit - 1);
        }
        vec2 p = 2*(gl_PointCoord - 0.5);
        if (markerShape == 1) // shape: .
        {
            float r = length(p);
            if (r > 1)
                discard;
#           ifndef BROKEN_GL_FRAG_COORD
            gl_FragDepth += projectionMatrix[3][2] * gl_FragCoord.w*gl_FragCoord.w
                            // TODO: Why is the factor of 0.5 required here?
                            * 0.5*modifiedPointRadius*sqrt(1-r*r);
#           endif
        }
        else if (markerShape == 2) // shape: o
        {
            float r = length(p);
            if (r > 1 || r < 1 - w)
                discard;
        }
        else if (markerShape == 3) // shape: x
        {
            w *= 0.5*sqrt2;
            if (abs(p.x + p.y) > w && abs(p.x - p.y) > w)
                discard;
        }
        else if (markerShape == 4) // shape: +
        {
            w *= 0.5;
            if (abs(p.x) > w && abs(p.y) > w)
                discard;
        }
    }
    fragColor = vec4(pointColor, 1);
}

#endif