```
library(habtools)
library(raster)
```

In this vignette, we give a brief overview of fractal dimension
methods in `habtools`

. There are currently five methods
summarized in the table below. Each method has inherent biases and so
should be used cautiously and consistently. For example, some methods
are sensitive to particular elements, like outliers, axis-scaling or
drop-offs, where others are invariant to them.

Method | Function | Data | Default `lvec` |
Notes |
---|---|---|---|---|

Height variation | `hvar` |
RasterLayer | 2-10 x resolution to extent | Invariant to height. Sensitive to non-fractal edges and outliers. |

Standard deviation | `sd` |
RasterLayer | 2-10 x resolution to extent | Invariant to height. Sensitive to non-fractal edges and outliers. |

Area | `area` |
RasterLayer or mesh3d | 2-4 x resolution to 1/8 extent | Sensitive to height; approaches 2 for surfaces with low height variation, relative to the extent. Less sensitive to non-fractal anomalies. |

Cube counting | `cubes` |
RasterLayer or mesh3d | 2-10 x resolution to extent | Sensitive to height; approaches 2 for for surfaces with low height variation, relative to the extent. More reliable for closed 3D object. Less sensitive to non-fractal anomalies. |

Box counting | `boxes` |
xy data.frame | 2-10 x resolution to extent | For 2D shapes only. Less sensitive to non-fractal anomalies. |

Generally, fractal dimension is estimated by measuring values (e.g., counts of cubes, heights, areas) across a sequence of different resolutions or length scales, and then finding the slope between these scales and measured values.

`lvec`

Selecting a sensible sequence of scales (i.e., the length vector
`lvec`

) is critical and should be the same when comparing
surfaces and objects. The `lvec`

should range between the
resolution and extent of the surface or mesh, or across a range of
scales relevant to the question being address. If you do not enter an
`lvec`

, each function will estimate one for you based on some
rules of thumb for that particular method (see Table). For example,
`hvar`

, `sd`

, `cubes`

and
`boxes`

will generate a sequence from 2- to 10-times the
resolution to the extent on a log2-scale. `area`

will
generate a sequence from 2 times the resolution to 1/8 the extent on the
log2-scale.

The `area`

method has increased bias with increased values
inside `lvec`

, because of the way surface area is calculated
at each scale (using Jenness, 2004). The method calculates surface area
of each cell by using the height of the cell and its 8 surrounding cells
to create a surface. For border cells, the cells outside of the DEM are
assumed to be of equal height as the border cells. This means that
border cells may have lower surface area estimates. If the ratio of
border cells to total cells inside a DEM is high, there may be an
underestimation of surface area, which is the case when values inside
`lvec`

are too close to the extent. For example, if
`lvec`

contains a scale of extent/2, all the cells are border
cells; if it contains extent/4, 75% of the cells are border cells. For
this reason, the `area`

method may be better suited for
larger areas, relative to the `lvec`

.

Below a fractal 128 x 128 terrain is simulated and fractal dimension
is estimated with the height variation method. Because no
`lvec`

was provided, one was generated as
`c(8, 16, 32, 64, 128)`

.

```
# simulate fractal terrain
<- sim_dem(L = 64, smoothness = 0.5)
surf
# fractal dimension using height variation method
fd(surf, method="hvar", plot = TRUE, diagnose = TRUE)
#> lvec is set to c(8, 16, 32, 64).
```

```
#> $D
#> [1] 2.478227
#>
#> $data
#> l h
#> 1 8 0.5603946
#> 2 16 0.8430506
#> 3 32 1.1709355
#> 4 64 1.6768663
#>
#> $lvec
#> [1] 8 16 32 64
#>
#> $D_vec
#> [1] 2.410824 2.526029 2.481894
#>
#> $var
#> [1] 0.05812521
#>
#> $method
#> [1] "hvar"
```

When `plot=TRUE`

, a plot of the surface is generated with
red rectangles representing `lvec`

to help visualize the
scales at which fractal dimension is estimated. When
`diagnostic=TRUE`

, the relationship between scales and mean
surface heights is plotted. The red line in the best-fit line for all
point and dashed lines are for each pair-wise scale with localized D
estimates. For the height range method, D is 3 minus the slope. The
diagnostic plot helps identify scales that are problematic (typically
the extremes of `lvec`

), which then have bias the estimate.
The `lvec`

is returned and can then be altered, entered
explicitly, and fractal dimension estimated again.

`hvar`

and `sd`

methods are not sensitive to
z-axis scaling, whereas `area`

and `cubes`

methods
are. Here are estimates for the relative flat surface we simulated with
`smoothness=0.5`

.

```
# height variation
fd(surf, method = "hvar")
#> lvec is set to c(8, 16, 32, 64).
#> [1] 2.478227
# standard deviation
fd(surf, method = "sd")
#> lvec is set to c(8, 16, 32, 64).
#> [1] 2.594688
# area
fd(surf, method = "area")
#> lvec is set to c(2, 4, 8).
#> [1] 2.001035
# cube counting
fd(surf, method = "cubes")
#> lvec is set to c(2, 4, 8, 16, 32, 64).
#> [1] 2.005364
```

`hvar`

and `sd`

estimates are in the vicinity
of D=2.5. However, `area`

and `cubes`

methods
result in estimates close to 2: a flat surface. If we transform the
z-scale to match the extent (L=128), we get the following.

```
<- surf
surf_z values(surf_z) <- values(surf_z) * ((64-1) / hr(surf_z))
hr(surf_z)
#> [1] 63
# height variation
fd(surf_z, method = "hvar")
#> lvec is set to c(8, 16, 32, 64).
#> [1] 2.478227
# standard deviation
fd(surf_z, method = "sd")
#> lvec is set to c(8, 16, 32, 64).
#> [1] 2.594688
# area
fd(surf_z, method = "area")
#> lvec is set to c(2, 4, 8).
#> [1] 2.361038
# cube counting
fd(surf_z, method = "cubes")
#> lvec is set to c(2, 4, 8, 16, 32, 64).
#> [1] 2.330863
```

`hvar`

and `sd`

estimates are identical (i.e.,
invariant to z-scaling) and `area`

and `cubes`

estimates are higher. This dichotomy raises an underlying philosophical
question about fractal surfaces: does the fractal dimension of an
surface change if it is transformed along a dimension?

The answer will likely depend on what you are analyzing. For example,
if you are examining landscapes, then `area`

and
`cubes`

methods will not be much use as they will tend to D=2
when surfaces a low relative to wide. Whereas, if you are examining
surfaces that are as high as they are wide, then the `area`

and `cubes`

methods might be a better choice. If you are
examining a closed 3D objects (i.e., watertight), then the
`cubes`

method is preferable as the 3D analogue to box
counting on closed shapes like the coastline of the UK.

```
# cube counting
fd(mcap, method = "cubes", plot=TRUE, diagnose=TRUE)
#> lvec is set to c(0.053, 0.106, 0.212, 0.423).
```

```
#> $D
#> [1] 2.315246
#>
#> $data
#> l n
#> 4 0.05291204 134
#> 3 0.10582408 31
#> 2 0.21164817 8
#> 1 0.42329634 1
#>
#> $lvec
#> [1] 0.42329634 0.21164817 0.10582408 0.05291204
#>
#> $D_vec
#> [1] 2.111893 1.954196 3.000000
#>
#> $var
#> [1] 0.5638126
#>
#> $method
#> [1] "cubes"
```

Here is an example of 2D box counting for the planar projection of
the `mcap`

mesh.

```
# project coral as xy coordinates
<- mesh_to_2d(mcap)
mcap_2d
# box counting
fd(mcap_2d, method = "boxes", plot=TRUE, diagnose=TRUE)
```

```
#> $D
#> [1] 1.39879
#>
#> $data
#> l n
#> 6 0.01219243 138
#> 5 0.02438486 65
#> 4 0.04876971 28
#> 3 0.09753943 11
#> 2 0.19507886 4
#> 1 0.39015772 1
#>
#> $lvec
#> [1] 0.39015772 0.19507886 0.09753943 0.04876971 0.02438486 0.01219243
#>
#> $D_vec
#> [1] 1.086157 1.215013 1.347923 1.459432 2.000000
#>
#> $var
#> [1] 0.3523516
#>
#> $method
#> [1] "boxes"
```

Natural surfaces often have anomalies such as peaks, troughs and
edges that can result in biased fractal dimension estimates. It is up to
the user to determine which method is best for their purposes, and to
demonstrate that a biological or ecological result is robust to the
choice of method. The function `detect_drop()`

can help to
investigate if a DEM has a lot of overhangs or sudden drops. It returns
a raster where values mark the areas with sudden drops depending on a
given threshold (Default threshold = 0.1).

`<- dem_crop(horseshoe, x0 = -470.8104, y0 = 1270.625, L = 2, plot = TRUE) dem1 `

```
<- detect_drop(dem1, d = 0.1)
drop1 plot(drop1)
```

```
# This DEM does not have many drops
fd(dem1, method = "hvar", lvec = c(1, 0.5, 0.25, 0.125), diagnose=TRUE)
```

```
#> $D
#> [1] 2.35037
#>
#> $data
#> l h
#> 1 0.125 0.09603534
#> 2 0.250 0.16825392
#> 3 0.500 0.24666635
#> 4 1.000 0.37923605
#>
#> $lvec
#> [1] 0.125 0.250 0.500 1.000
#>
#> $D_vec
#> [1] 2.190997 2.448079 2.379465
#>
#> $var
#> [1] 0.1331159
#>
#> $method
#> [1] "hvar"
dem2 <- dem_crop(horseshoe, x0 = -466.8104, y0 = 1266.625, L = 2, plot = TRUE)
```

```
<- detect_drop(dem2, d = 0.1)
drop2 plot(drop2)
```

`fd(dem2, method = "hvar", lvec = c(1, 0.5, 0.25, 0.125), diagnose=TRUE)`

```
#> $D
#> [1] 2.198017
#>
#> $data
#> l h
#> 1 0.125 0.1849993
#> 2 0.250 0.3998837
#> 3 0.500 0.7186458
#> 4 1.000 0.9706139
#>
#> $lvec
#> [1] 0.125 0.250 0.500 1.000
#>
#> $D_vec
#> [1] 1.887939 2.154299 2.566383
#>
#> $var
#> [1] 0.3418205
#>
#> $method
#> [1] "hvar"
```

In the example above, `dem2`

has much more and higher
drops compared to `dem1`

. We can see its effect on the height
variation method by looking at the diagnostics plot. The dots in the
diagnostic plot for `dem2`

follow a curve, and do not fall on
a straight line. The variance of D in `dem2`

(var = 0.44) is
also larger than for `dem1`

(0.24). Moreover, the frequency
of drops correlates slightly negatively with the fractal dimension
estimates when using the `hvar`

method. Therefore, caution
must be used with the interpretation of fractal dimension when the DEM
has many drops or overhangs. Changes in D depending on the scale is
common for non-fractal surfaces and therefore it is crucial to carefully
choose the `lvec`

and keep it the same when comparing across
many surfaces.

- Jenness, J.S. Calculating Landscape Surface Area from Digital Elevation Models. Wildlife Society Bulletin, Vol. 32, No. 3 (Autumn, 2004), pp. 829-839.