Get geodesic distortion for rasters

Calculate the geodesic distortion on a raster

getDistortion(z)

Arguments

z

An input SpatRaster with a known projection

Value

A SpatRaster with three layers: (1) 'lx' with the horizontal pixel size

(2) 'ly' with the vertical pixel size

(3) 'A' with the pixel area

Examples

# Generate dummy dem, assign it a projection
n <- 5
dem <- expand.grid(list(x = 1:(n * 100),
                        y = 1:(n * 100))) / 100
dem <- as.data.table(dem)
dem[, z := 250 * exp(-(x - n/2)^2) + 
      250 * exp(-(y - n/2)^2)]
#>            x    y         z
#>      1: 0.01 0.01 1.0146139
#>      2: 0.02 0.01 1.0404641
#>      3: 0.03 0.01 1.0675194
#>      4: 0.04 0.01 1.0958300
#>      5: 0.05 0.01 1.1254477
#>     ---                    
#> 249996: 4.96 5.00 1.0711366
#> 249997: 4.97 5.00 1.0428260
#> 249998: 4.98 5.00 1.0157707
#> 249999: 4.99 5.00 0.9899205
#> 250000: 5.00 5.00 0.9652271
dem <- rast(dem)
ext(dem) <- c(10000, 20000, 30000, 40000)
crs(dem) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84"

# Get the geodesic distorsion
distort <- getDistortion(dem)