-
Notifications
You must be signed in to change notification settings - Fork 4
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- BPI Added - WoodEvans now calculates the representative measures of plan, profile, and twisting curvature (Minar et al., 2020) - WoodEvans crosc and longc removed - - WoodEvans replaced directional derivative versions of minc, maxc, and meanc with kmin, kmean, kmax (Minar et al., 2020) - WoodEvans features algorithm modified to use kmin, kmax, and planc
- Loading branch information
Showing
14 changed files
with
369 additions
and
67 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,130 @@ | ||
| #' Creates circular focal window | ||
| #' | ||
| #' Creates circular focal window around central pixel. | ||
| #' @param radius radius of circular window | ||
| #' @param resolution resolution of intended raster layer (one number or a vector of length 2). Only necessary if unit= "map" | ||
| #' @param unit unit for radius. Either "cell" (number of cells, the default) or "map" for map units (e.g. meters). | ||
| #' @param return_dismat logical, if TRUE return a matrix of distances from focal cell instead of a matrix to pass to terra::focal | ||
| #' @export | ||
| circle_window<- function(radius, unit= "cell", resolution, return_dismat = FALSE){ | ||
| if (unit != "map" & unit != "cell"){ | ||
| stop("Error: unit must equal 'map' or 'cell'") | ||
| } | ||
| if(length(radius)!=1){ | ||
| stop("Error: radius must be a single integer") | ||
| } | ||
| if(unit=="cell"){ | ||
| resolution<- c(1,1) | ||
| } | ||
| if (length(resolution)==1){ | ||
| resolution<- rep(resolution,2) | ||
| } | ||
| nrows<- floor((radius/resolution[2])*2+1) | ||
| ncols<- floor((radius/resolution[1])*2+1) | ||
| if(nrows %% 2 ==0){ | ||
| nrows<- nrows+1 | ||
| } | ||
| if(ncols %% 2 ==0){ | ||
| ncols<- ncols+1 | ||
| } #nrow and ncol must be odd to have central pixel | ||
| x<- matrix(seq(1:ncols), nrow = nrows, ncol =ncols, byrow = TRUE) - ((ncols+1)/2) | ||
| y<- matrix(seq(1:nrows), nrow=nrows, ncol=ncols, byrow = FALSE) - ((nrows+1)/2) | ||
| x<- x * resolution[1] | ||
| y<- y * resolution[2] | ||
| dis_mat<- sqrt((y^2)+(x^2)) #Distance from center of window | ||
| if(return_dismat){ | ||
| return(dis_mat) | ||
| } | ||
| w<- matrix(NA, nrow = nrow(dis_mat), ncol= ncol(dis_mat)) | ||
| w[dis_mat <= radius]<- 1 | ||
| return(w) | ||
| } | ||
|
|
||
| #' Creates annulus focal window | ||
| #' | ||
| #' Creates annulus focal window around central pixel. | ||
| #' @param radius radius of inner annulus c(inner,outer) | ||
| #' @param unit unit for radius. Either "cell" (number of cells, the default) or "map" for map units (e.g. meters). | ||
| #' @param resolution resolution of intended raster layer (one number or a vector of length 2). Only necessary if unit= "map" | ||
| #' @param return_dismat logical, if TRUE return a matrix of distances from focal cell instead of a matrix to pass to terra::focal (default FALSE) | ||
| #' @export | ||
| annulus_window<- function(radius, unit= "cell", resolution, return_dismat=FALSE){ | ||
| if(length(radius)==1){radius<- rep(radius, 2)} | ||
| if(length(radius) > 2){ | ||
| stop("Specified radius exceeds 2 dimensions") | ||
| } | ||
| if(radius[1] > radius[2]){ | ||
| stop("Error: inner radius must be less than or equal to outer radius") | ||
| } | ||
| if((radius[1] < 1) & (unit == "cell")){ | ||
| stop("Error: inner radius must be at least 1") | ||
| } | ||
| if(unit=="cell"){ | ||
| resolution<- c(1,1) | ||
| } | ||
| if((radius[1] < max(resolution)) & (unit == "map")){ | ||
| stop("Error: inner radius must be >= resolution") | ||
| } | ||
| if(length(resolution) > 2){ | ||
| stop("Specified inner radius exceeds 2 dimensions") | ||
| } | ||
|
|
||
| dis_mat<- circle_window(radius = radius[2], unit = unit, resolution = resolution, return_dismat = TRUE) | ||
| if(return_dismat){ | ||
| return(dis_mat) | ||
| } | ||
| w<- matrix(NA, nrow=nrow(dis_mat), ncol=ncol(dis_mat)) | ||
| w[(dis_mat >= radius[1]) & (dis_mat <= radius[2])]<- 1 | ||
| return(w) | ||
| } | ||
|
|
||
| #' Calculates Bathymetric Position Index | ||
| #' | ||
| #' Calculates Bathymetric Position Index (BPI). This is the value of the focal pixel minus the mean of the surrounding pixels contained within an annulus shaped window. | ||
| #' @param r DEM as a SpatRaster or RasterLayer | ||
| #' @param radius a vector of length 2 specifying the inner and outer radii of the annulus c(inner,outer). This is ignored if w is provided. | ||
| #' @param unit unit for radius. Either "cell" (number of cells, the default) or "map" for map units (e.g. meters). This is ignored if w is provided. | ||
| #' @param w A focal weights matrix representing the annulus focal window created using MultiscaleDEM::annulus_window. | ||
| #' @param na.rm A logical vector indicating whether or not to remove NA values before calculations | ||
| #' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE). Only valid if radius is used. | ||
| #' @return a SpatRaster or RasterLayer | ||
| #' @import terra | ||
| #' @importFrom raster raster | ||
| #' @references | ||
| #' Lundblad, E.R., Wright, D.J., Miller, J., Larkin, E.M., Rinehart, R., Naar, D.F., Donahue, B.T., Anderson, S.M., Battista, T., 2006. A benthic terrain classification scheme for American Samoa. Marine Geodesy 29, 89–111. | ||
| #' @export | ||
|
|
||
| BPI<- function(r, radius=NULL, unit= "cell", w=NULL, na.rm=FALSE, include_scale=FALSE){ | ||
| og_class<- class(r)[1] | ||
| if(og_class=="RasterLayer"){ | ||
| r<- terra::rast(r) #Convert to SpatRaster | ||
| } | ||
| #Input checks | ||
| if(!(og_class %in% c("RasterLayer", "SpatRaster"))){ | ||
| stop("Error: Input must be a 'SpatRaster' or 'RasterLayer'") | ||
| } | ||
| if(terra::nlyr(r)!=1){ | ||
| stop("Error: Input raster must be one layer.") | ||
| } | ||
| if(!is.null(w)){ | ||
| if(class(w)[1] != "matrix"){ | ||
| stop("Error: w must be a matrix") | ||
| } | ||
| if(any(0 == (dim(w) %% 2))){ | ||
| stop("Error: dimensions of w must be odd") | ||
| } | ||
| } | ||
| if(length(radius)==1){radius<- rep(radius, 2)} | ||
|
|
||
| resolution<- terra::res(r) | ||
|
|
||
| if(is.null(w)){ | ||
| w<- annulus_window(radius = radius, unit = unit, resolution = resolution) | ||
| } | ||
|
|
||
| bpi<- r - terra::focal(x = r, w = w, fun = mean, na.rm = na.rm) | ||
| names(bpi)<- "BPI" | ||
| if(include_scale & (!is.null(radius))){names(bpi)<- paste0(names(bpi), "_", radius[1],"x", radius[2])} #Add scale to layer names | ||
| if(og_class=="RasterLayer"){bpi<- raster::raster(bpi)} | ||
| return(bpi) | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,41 @@ | ||
| knc<- function(a,b,c,d,e,f){ | ||
| #Planc: Normal Contour Curvature | ||
| out<- -(2*a*(e^2) - 2*c*d*e + 2*b*d^2)/((d^2+e^2) * sqrt(1+d^2+e^2)) | ||
| return(out) | ||
| } | ||
|
|
||
| kns<- function(a,b,c,d,e,f){ | ||
| #Profc: Normal Slope Line Curvature | ||
| out<- (-2 * (a*d^2 + b*e^2 + c*d*e)) / ((e^2 + d^2)*(1 + e^2 + d^2)^1.5) | ||
| return(out) | ||
| } | ||
|
|
||
| tgc<- function(a,b,c,d,e,f){ | ||
| #TwistC: Contour geodesic torsion | ||
| out<- (d*e*(2*a-2*b) - c*(d^2-e^2))/((d^2+e^2)*(1+d^2+e^2)) | ||
| return(out) | ||
| } | ||
|
|
||
| kmean<- function(a,b,c,d,e,f){ | ||
| #Mean Curvature | ||
| out<- -((1+e^2)*2*a - 2*c*d*e + (1+d^2)*2*b) / (2*sqrt((1+d^2+e^2)^3)) | ||
| return(out) | ||
| } | ||
|
|
||
| ku<- function(a,b,c,d,e,f){ | ||
| #unsphericity curvature | ||
| out<- sqrt((((1+e^2)*2*a - 2*c*d*e +(1+d^2)*2*b) / (2*sqrt((1+d^2+e^2)^3)))^2 - ((2*a*2*b-c^2)/(1+d^2+e^2)^2)) | ||
| return(out) | ||
| } | ||
|
|
||
| kmin<- function(a,b,c,d,e,f){ | ||
| #Min Curvature | ||
| out<- kmean(a,b,c,d,e)-ku(a,b,c,d,e) | ||
| return(out) | ||
| } | ||
|
|
||
| kmax<- function(a,b,c,d,e,f){ | ||
| #Max Curvature | ||
| out<- kmean(a,b,c,d,e)+ku(a,b,c,d,e) | ||
| return(out) | ||
| } |
Oops, something went wrong.