ailich / multiscaledtm Goto Github PK

Update code from using raster to terra

Many raster functions have a filename argument where you can write a raster to file and also have it available in the R environment. Add this feature to functions.

add examples to documentation

According to Wilson et al 2007, y should decrease as you go down rows (not increase). See how this affects calculated metrics as this affects the sign of regression parameters c and e.

"the use of the terms DEM, DBM, CTM, and DDM should be abandoned in favour of the more neutral, all-encompassing term DTM" (Lecours et al., 2016)

Options to pass to terra. WIll make debugging easier.

Will be easier with terra which supports NA focal weights

"Error: cannot allocate vector of size" or "Error: memory exhausted (limit reached?)" when working on large rasters.

library(terra)
library(MultiscaleDEM)
library(GLCMTextures)

stk_full<- rast("Data/namibia_full.tif")
grn <- stk_full$green

grn_quantized <- quantize_raster(r = grn, n_levels = 16,
                                 method = "equal prob")

out<- Qfit(grn_quantized, w=c(5,5), metrics = c(), na.rm=TRUE, return_params = TRUE)

See ailich/GLCMTextures#12 and rspatial/terra#479

9

Add option to choose which metrics are calculated. Maybe even with negative indexing (e.g. specifying which metrics not to calculate). Have it so metrics are only calculated if needed, not all calculated and then subset.

In certain cases getting error related to matrix being singular (can't be inverted)

#bathy is BTM raster
SDAdj<- AdjSD(bathy, w=c(7,7), na.rm = TRUE, pad = TRUE)
#error: inv(): matrix is singular

Add cff citation file (https://citation-file-format.github.io/)

After doing some checking and comparing with Minar et al 2020, Wood 1996 and r.param.scale have a sign convention for plan curvature that is inconsistent with the other curvatures. Flip sign convention

AdjSD should use x and y but not xy. xy allows for curvature and therefore doesn't fit a plane to the data.

"The grid of residual topography was then calculated as the cell by-cell difference between the LiDAR DTM and the mean DTM. Finally, the standard deviations of residual topography values are computed" (Cavalli et al 2008)

Remove intercept column of ones from X matrix and subtract elevation value of focal cell so that the focal cell is located at 0,0,0 in the local coordinate system

xy resolution is flipped.
PA<- (x_res*w[1]) * (y_res*w[2]) should be PA<- (x_res*w[2]) * (y_res*w[1]) because width is x_resnumber of columns and height is y_resnumber of rows.

Add option to set extreme outliers in Qfit parameters to NA

Losing factor levels in features when converting to raster in WoodEvans

Add Terrain Attributes Explorer shiny app to package where it can be run using a function.

https://ailich.shinyapps.io/Terrain_Attributes_Explorer_App/

https://github.com/ailich/Terrain_Attributes_Explorer

https://deanattali.com/2015/04/21/r-package-shiny-app/

Add na.rm option to VRM function. Adjust focal(..., fun=sum) to have na.rm in it. When na.rm=TRUE, adjust scale.factor (n) to count up only the number of cells that aren't NA.

Determine and implement a list of test cases to make sure functions are working as intended after each update

Getting stripes between chunks when merging rasters too large to fit in memory for AdjSD and WoodEvans

Document which functions can only use projected coordinate systems and throw an error if an unprojected coordinate system is used (i.e. lat/lon).

Use raster::isLonLat/raster::couldBeLonLat

Remove arma::mat H = X * XtX_inv * Xt from C_OLS_params. Don't need H matrix to calculate regression coefficients.

I can't find the full references anywhere for the sources you cite in various functions. I think the documentation for each function should contain the relevant references.

Add TPI and RDMV

The formulas for these two are the same in Wood 1996 pages 86 (prof curv) and page 115 (max/min curv). This may be describing only a special case. max/min curv is described as the maximum/minimum curvature (in any plane) so maybe use calc across a stack of curvatures and apply max and min as the functions.

Qfit currently will mask values if any regression parameter =0 rather than if all are equal to zero. Better to do sum(abs(params))==0

Should the scale of the slope layer be w or w+2?

• Add option to output derivatives (zxx, zyy, zxy, zx, zy) instead of regression coefficients (a-f) so that Minar et al 2020 formulas can be used directly.
• Use the basic trio of curvatures: profile, plan and twisting curvatures.
  • Add a measure of "twisting curvature"
  • Use the "principal representative" for each measure of curvature in the basic trio
    • "The principal representatives (kn)s [profile/normal slope line curvature], (kn)c [plan/normal contour curvature] and (τg)c [twisting/contour geodesic torsion] have demonstrably greater interpretable robustness than other curvatures of their particular type." (Minar et al 2020)
• Change mean curvature from z''mean (Mean 2nd derivative) to kmean
  • "Wilson et al. (2007) used these zmax′ ′ and zmin′′ for computation of a ‘mean curvature’ that is in fact the mean second directional derivative of altitude (zmean′ ′) – double the Elevation Laplacian (∇2z) and significantly different from the geometric kmean." (Minar et al 2020)

Add option to fit a 3rd order instead of second order polynomial in Qfit (Florinsky, 2016).

Florinsky, I.V., 2016. Digital terrain analysis in soil science and geology, 2nd ed. Elsevier/Academic Press, Amsterdam.

Remove crosc and longc from default metrics as they are no longer included as one of the available metrics

"The regular nature of the DEM grid requires only one inversion of the matrix for any combination of window size and weight matrix. This is computationally very advantageous in the application of quadratic surface modelling over large DEMs."- Albani et al., 2004

However, code currently uses Armadillo's solve function b/c of numeric instabilities related to matrix inversion.

Maybe split Planar feature into "FLAT" (slope <= slope threshold) and "SLOPE" (slope > slope threshold)

Have option to have weights for each cell (e.g. inverse distance weighted).

From Wood, 1996:
"In common with all inverse distance functions, an (arbitrary) decision has to be made about the weight of the central cell, which by strict definition would have an undefined weight of 1/0. To avoid the problem of infinite weights, unity is added to each distance in the following weighting function: wij = 1/(dij+1)^n where dy is the Euclidean distance in grid cells to the central cell and n is an exponent ranging from 0 (no distance decay), through 1 (linear decay), to 2 (distance squared decay)."

Plan curvature in WoodEvans tends to have a few very extreme values (particularly at smaller window sizes). Maybe a TRUE/FALSE toggle specifying to remove "extreme values" or a numeric representing a user specified proportion of the data to remove at both ends which can be supplied to raster::quantile

thresh<- 0.01 #remove upper and lower 1% of data values
min_thresh<- raster::quantile(plancurv, thresh)
max_thresh<- raster::quantile(plancurv, 1-thresh)
plancurv[plancurv< min_thresh]<- NA
plancurv[plancurv > max_thresh]<- NA

Error in quant[, 3] : incorrect number of dimensions. Related to subsetting a matrix vs vector in outlier detection code.

in cpp_code out mask layer should be initialized as out(_,6)=rep(0,out.nrow()) not out(_,6)=rep(0,100) but it does not affect results as the whole column still fills in with 0's.

To remove correlation among successively larger spatial scales that include data from all previous scales, add option for the quadratic to be fit only using only the outer edge of the focal window and the central cell (similar to what SlpAsp does).

devtools::check results

-- R CMD check results ---------------------------------------------------------------------------------------------------- MultiscaleDEM 0.0.2.1 ----
Duration: 1m 53.4s

> checking Rd \usage sections ... WARNING
  Undocumented arguments in documentation object 'SAPA'
    'slope_correction'
  
  Functions with \usage entries need to have the appropriate \alias
  entries, and all their arguments documented.
  The \usage entries must correspond to syntactically valid R code.
  See chapter 'Writing R documentation files' in the 'Writing R
  Extensions' manual.

> checking DESCRIPTION meta-information ... NOTE
  Malformed Description field: should contain one or more complete sentences.
  Non-standard license specification:
    `use_mit_license()`
  Standardizable: FALSE

> checking top-level files ... NOTE
  File
    LICENSE
  is not mentioned in the DESCRIPTION file.
  Non-standard files/directories found at top level:
    'README.Rmd' 'README_files' 'images'

0 errors √ | 1 warning x | 2 notes x

Make code more memory safe by using functions instead of indexing by a raster layer (rspatial/terra#325)

Seems like weights can be used to account for if there are missing values in the window.

From ArcGIS How Slope Works
" [dz/dx] = ((c + 2f + i)*4/wght1 - (a + 2d + g)*4/wght2) / (8 * x_cellsize)
where wght1 and wght2 are the horizontal weighted counts of valid cells.

For instance, if c, f, and i all have valid values, wght1 = (1+21+1) = 4.
i is NoData, wght1 = (1+21+0) = 3.
f is NoData, wght1 = (1+2*0+1) = 2.

Similar logic applies to wght2, except the neighbor locations are a, d, and g.

The rate of change in the y direction for cell e is calculated with the following algorithm:

[dz/dy] = ((g + 2h + i)*4/wght3 - (a + 2b + c)*4/wght4) / (8 * y_cellsize)

where wght3 and wght4 are the same concept as in the [dz/dx] computation."

Related to issue in raster package rspatial/raster#228

Explore adding support for more functions to be used in unprojected coordinate systems.

raster::terrain (terrain.R, terrain.cpp) uses haversine distance and ArcPro Surface Parameters uses a 3D Earth Centered, Earth Fixed coordinate system.