AustralianAntarcticDivision / rema.proc

processing Antarctic Reference Elevation tiles (REMA) at various resolutions (8m, 100m, 200m)
3 stars 0 forks source link

Compass-based aspect #5

Open raymondben opened 3 years ago

raymondben commented 3 years ago

This will convert aspect-relative-to-grid (i.e. the current aspect tifs) to aspect-relative-to-compass, which is probably what most users will want:

aspect_to_compass <- function(x) {
    crds <- coordinates(x)
    ## calculate longitude of each grid cell in radians
    gridlon <- pi/2 - atan2(crds[, 2], crds[, 1])
    angle_normalise <- function(x) (x + pi) %% (2 * pi) - pi ## normalize angle in radians to range [-pi,pi)
    temp <- angle_normalise(values(x) / 180 * pi) ## convert aspect to radians and normalise to -pi,pi
    temp <- temp - gridlon ## adjust aspect for longitude
    temp <- angle_normalise(temp) / pi * 180 ## normalize again and back to degrees
    values(x) <- temp
    x
}

I've run it for the 1km and 200m single-tif files (see *_aspect_compass.tif). Probably worth pre-computing at 100m as well, but then bundling this little helper into raadtools for users to apply on the fly to 8m? (Would be wise to add a check to it, to ensure that the input x is in zero-centred polar stereo, otherwise it willna work)

mdsumner commented 3 years ago

look at this

https://gis.stackexchange.com/a/203109/1204

...