This is building off of the first example from Sierra’s markdown that she shared via email on 1/24/2024. It just so happens that it’s my trajectory :).

library(lidR)
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
library(terra)
## terra 1.7.46
# read in the data
las <- readLAS("S:/ursa/cutler/transects/Hovermap_copies/singlepoints/T5_MJC_traj.las", select = "xyz")
dtm <- rast("S:/ursa2/campbell/hvm_travel_rate/data/als/data/dtm_01m.tif")

# extract dtm values at trajectory points
las <- merge_spatial(las, dtm, "Z.dtm")

# adjust for height
las$Z <- las$Z - 1.53

# get difference
las <- add_lasattribute(las, las$Z - las$Z.dtm, "Z.diff", "HVM elev minus DTM elev")

# plot differences in x-y space
n <- 100
pal <- colorRampPalette(c("red", "lightgray", "blue"))(n)
max.val <- max(abs(las$Z.diff))
min.val <- max.val * -1
brks <- seq(min.val, max.val, length.out = n)
cols <- pal[cut(las$Z.diff, breaks = brks)]
plot(las$X, las$Y, pch = 16, col = cols)
legend("topleft", legend = c(round(min.val,2), 0, paste0("+", round(max.val,2))),
       pch = 16, col = c("red", "lightgray", "blue"))

Cool, so it looks like HVM elevations are underestimated on the SW portion of the trajectory, and overestimated on the NE portion. We want to be able to correct for those local differences when adjusting HVM elevations to match USGS elevations. That’s, of course, easy enough to do on the trajectory, but we need to do this not just on but also in the area surrounding the trajectory. So we’ll need to interpolate a surface that represents a per-pixel correction factor that should enable us to correct all points in the HVM data. There are many different ways we could go about this, but I think a key point is that we don’t want the interpolation to be too detailed. Clearly the trajectories feature the undulation of up and down stepping motions, so if we fit a super detailed correction surface, it would include those undulations, which we obviously don’t want. We just want to capture the spatial trend in differences between elevations.

I think a simple way to start would be a linear or polynomial model.

# define variables of interest
x <- las$X
y <- las$Y
z <- las$Z.diff

# fit linear 2D model (1st order polynomial)
mod.lm.1 <- lm(z ~ x + y)
summary(mod.lm.1)
## 
## Call:
## lm(formula = z ~ x + y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44619 -0.08369  0.00870  0.08214  0.42706 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.691e+04  6.466e+02 103.468  < 2e-16 ***
## x           -2.126e-03  3.308e-04  -6.428  1.3e-10 ***
## y           -1.467e-02  1.774e-04 -82.702  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1238 on 76195 degrees of freedom
## Multiple R-squared:  0.9382, Adjusted R-squared:  0.9382 
## F-statistic: 5.786e+05 on 2 and 76195 DF,  p-value: < 2.2e-16

…Pretty decent performance (94% of variance explained). Let’s map it out using terra::interpolate(). But, to do that, we need to supply a reference raster, whose extent, coordinate system, pixel grid, etc. will be matched to. We don’t want to use the entire DTM, since it’s massive. So, we’ll just buffer our the trajectory by 10m, clip the DTM to that, and use that as the reference first.

# convert xy points to SpatVector
pts <- vect(las@data, crs = wkt(crs(las)), geom = c("X", "Y"))

# buffer them
buff <- buffer(pts, 10) |> aggregate()

# clip the dtm
dtm.clip <- crop(dtm, buff, mask = T)

# plot everything out
plot(dtm.clip)
lines(buff, col = "red", lwd = 2)
points(pts, col = "blue")

Right-o! Now back to interpolation…

# interpolate 1st order polynomial
r.lm.1 <- interpolate(dtm.clip, mod.lm.1) |> crop(buff, mask = T)

# plot it out
plot(r.lm.1)

Looks… Pretty good? Tough to say. What we can do to test it would be:

  1. Extract the values of this “correction” raster to the trajectory
  2. Subtract those values from the HVM elevations
  3. Compare corrected HVM elevations to USGS elevations to see how well they line up

LET’S GOOOOO

# extract interpolated correction raster values for trajectory points
las <- merge_spatial(las, r.lm.1, "Z.corr")

# add correction factor to HVM values
las <- add_lasattribute(las, las$Z - las$Z.corr, "Z.adj", "Adjusted elevation value")

# check out the first few rows
head(las@data)
##           X       Y        Z    Z.dtm   Z.diff   Z.corr    Z.adj
## 1: 457143.5 4493341 2532.872 2531.573 1.299342 1.218061 2531.654
## 2: 457143.5 4493341 2532.872 2531.573 1.299342 1.218061 2531.654
## 3: 457143.5 4493341 2532.872 2531.573 1.299342 1.218061 2531.654
## 4: 457143.5 4493341 2532.872 2531.573 1.299342 1.218061 2531.654
## 5: 457143.5 4493341 2532.872 2531.573 1.299342 1.218061 2531.654
## 6: 457143.5 4493341 2532.872 2531.573 1.299102 1.218061 2531.654

In the table above, Z is the raw HVM elevations, Z.dtm is the USGS elevations, and Z.adj is the new, interpolation-surface-adjusted values. Pretty good start. Let’s take a look at a few histograms and some numbers:

# calculate new difference between Z.adj and Z.dtm
las <- add_lasattribute(las, las$Z.adj - las$Z.dtm, "Z.adj.diff", "Adjusted elevation value diff")

# generate histogram of differences
hist(las$Z.adj.diff)
mn <- mean(las$Z.adj.diff)
sd <- sd(las$Z.adj.diff)
abline(v = mn, lwd = 2, col = "blue")
abline(v = mn + sd, lwd = 2, col  = "red")
abline(v = mn - sd, lwd = 2, col  = "red")
legend("topleft", legend = c(paste0("mean = ", round(mn, 2)), paste0("sd = ", round(sd, 2))), lwd = 2, col = c("blue", "red"))

Pretty good! Let’s see how much better a second order polynomial does…

# fit 2nd order polynomial
mod.lm.2 <- lm(z ~ poly(x,2) + poly(y,2))
summary(mod.lm.2)
## 
## Call:
## lm(formula = z ~ poly(x, 2) + poly(y, 2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44785 -0.08373  0.00903  0.07744  0.42496 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  5.012e-01  4.482e-04 1118.080  < 2e-16 ***
## poly(x, 2)1 -1.180e+01  1.609e+00   -7.329 2.34e-13 ***
## poly(x, 2)2 -6.696e-01  6.930e-01   -0.966  0.33393    
## poly(y, 2)1 -1.215e+02  1.613e+00  -75.306  < 2e-16 ***
## poly(y, 2)2  2.104e+00  6.848e-01    3.072  0.00212 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1237 on 76193 degrees of freedom
## Multiple R-squared:  0.9383, Adjusted R-squared:  0.9383 
## F-statistic: 2.898e+05 on 4 and 76193 DF,  p-value: < 2.2e-16
# interpolate 2nd order polynomial
r.lm.2 <- interpolate(dtm.clip, mod.lm.2) |> crop(buff, mask = T)

# plot it out
plot(r.lm.2)

# extract interpolated correction raster values for trajectory points
las <- merge_spatial(las, r.lm.2, "Z.corr")

# add correction factor to HVM values
las <- add_lasattribute(las, las$Z - las$Z.corr, "Z.adj", "Adjusted elevation value")

# check out the first few rows
head(las@data)
##           X       Y        Z    Z.dtm   Z.diff  Z.corr    Z.adj Z.adj.diff
## 1: 457143.5 4493341 2532.872 2531.573 1.299342 1.22543 2531.647 0.08128075
## 2: 457143.5 4493341 2532.872 2531.573 1.299342 1.22543 2531.647 0.08128075
## 3: 457143.5 4493341 2532.872 2531.573 1.299342 1.22543 2531.647 0.08128075
## 4: 457143.5 4493341 2532.872 2531.573 1.299342 1.22543 2531.647 0.08128075
## 5: 457143.5 4493341 2532.872 2531.573 1.299342 1.22543 2531.647 0.08128075
## 6: 457143.5 4493341 2532.872 2531.573 1.299102 1.22543 2531.647 0.08104075
# calculate new difference between Z.adj and Z.dtm
las <- add_lasattribute(las, las$Z.adj - las$Z.dtm, "Z.adj.diff", "Adjusted elevation value diff")

# generate histogram of differences
hist(las$Z.adj.diff)
mn <- mean(las$Z.adj.diff)
sd <- sd(las$Z.adj.diff)
abline(v = mn, lwd = 2, col = "blue")
abline(v = mn + sd, lwd = 2, col  = "red")
abline(v = mn - sd, lwd = 2, col  = "red")
legend("topleft", legend = c(paste0("mean = ", round(mn, 2)), paste0("sd = ", round(sd, 2))), lwd = 2, col = c("blue", "red"))

Nearly identical performance to the first order. So, at least in this very specific case of this very specific trajectory, a 1st order polynomial interpolation correction does a good job of adjusting HVM elevations to match USGS elevations. The degree to which that will be the case on other trajectories remains to be seen. The spatial trend in differences between the two elevation measures was clearly very linear in this case. If that trend is not found on other trajectories, then perhaps a higher-order polynomial function would be better suited. Or, if the pattern is really spatially complex, we could even explore more advanced interpolation techniques such as inverse distance weighting or kriging. But, my gut says those might be overkill for this type of scenario.

Moving forward with the 1st order poly approach, here’s a quick comparison of the pre- and post-correction elevation differences, first in histogram form:

# refit 1st order polynomial
mod.lm.1 <- lm(z ~ x + y)

# interpolate 2nd order polynomial
r.lm.1 <- interpolate(dtm.clip, mod.lm.1) |> crop(buff, mask = T)

# extract interpolated correction raster values for trajectory points
las <- merge_spatial(las, r.lm.1, "Z.corr")

# add correction factor to HVM values
las <- add_lasattribute(las, las$Z - las$Z.corr, "Z.adj", "Adjusted elevation value")

# calculate new difference between Z.adj and Z.dtm
las <- add_lasattribute(las, las$Z.adj - las$Z.dtm, "Z.adj.diff", "Adjusted elevation value diff")

# plot combined histogram
h.raw <- hist(las$Z.diff, plot = F)
h.cor <- hist(las$Z.adj.diff, plot = F)
xrange <- range(c(las$Z.diff, las$Z.adj.diff))
yrange <- range(c(h.raw$counts, h.cor$counts))
plot(h.raw, main = NA, xlab = "Z Difference", border = NA, col = rgb(1,0,0,0.5),
     xlim = xrange, ylim = yrange)
plot(h.cor, main = NA, xlab = NA, border = NA, col = rgb(0,0,1,0.5),
     xlim = xrange, ylim = yrange, add = T)
legend("topright", legend = c("Uncorrected", "Corrected"),
       fill = c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)))

…Obviously much better. Now let’s take a look at the two mapped out to see how the differences play out in x-y space. The goal is that we no longer see a prevailing spatial trend in differences.

# plot differences in x-y space
cols.raw <- pal[cut(las$Z.diff, breaks = brks)]
cols.cor <- pal[cut(las$Z.adj.diff, breaks = brks)]
par(mfrow = c(1,2))
plot(las$X, las$Y, pch = 16, col = cols.raw, main = "Uncorrected")
legend("topleft", legend = c(round(min.val,2), 0, paste0("+", round(max.val,2))),
       pch = 16, col = c("red", "lightgray", "blue"))
plot(las$X, las$Y, pch = 16, col = cols.cor, main = "Corrected")
legend("topleft", legend = c(round(min.val,2), 0, paste0("+", round(max.val,2))),
       pch = 16, col = c("red", "lightgray", "blue"))

A clear improvement. There do appear to be some clusters of red in there, so I’m going to re-plot that one with a color ramp range that more clearly highlights contrast in the corrected values.

# plot differences in x-y space
max.val <- max(abs(las$Z.adj.diff))
min.val <- max.val * -1
brks <- seq(min.val, max.val, length.out = n)
cols.cor <- pal[cut(las$Z.adj.diff, breaks = brks)]
plot(las$X, las$Y, pch = 16, col = cols.cor, main = "Corrected")
legend("topleft", legend = c(round(min.val,2), 0, paste0("+", round(max.val,2))),
       pch = 16, col = c("red", "lightgray", "blue"))

A few interesting little patches of over- and underestimation. Not sure what to make of them. Again, it’s worth noting that these include the vertical undulations of stepping, climbing over and ducking under things. So, a few decimeters one direction or another may not be that problematic. But, if we do think this is a problematic level of over/underestimation, then we could explore a more spatially explicit interpolation technique, such as IDW.