README.md

area

The goal of area is to calculate areas, allow control over how that happens, and illustrate how it works.

Very much work in progress.

Please, do not use this for your land surveying contract.

Installation

You can install the development version from GitHub with:

# install.packages("devtools")
devtools::install_github("hypertidy/area")

Straightforward tools

The area package is written to allow easy to use calculations for developing packages. It’s similar to the traipse package which was created to leverage the grouping mechanisms of the tidyverse for common tracking data calculations.

library(area)
polygon_area(cbind(c(0, 1, 1, 0, 0), 
                c(0, 0, 1, 1, 0)))
#> [1] 1

triangle_area(cbind(c(0, 1, 1, 1, 0, 0), 
               c(0, 0, 1, 1, 1, 0)))
#> [1] 0.5 0.5

There’s an inbuilt data set with a path of Tasmania’s outline.

index <- 4423 ## explain later
idx <- c(tas$S[1], tas$S[1:index,2], 1)
plot(tas$P[idx, ], type = "l", asp = 1)

polygon_area(tas$P[idx, ])/1e6
#> [1] 61672.74

sf::st_area(sf::st_polygon(list(tas$P[idx, ])))/1e6
#> [1] 61672.74

library(dplyr)
library(area)

calc_area <- function(df) {
 area <- df %>% mutate(row = row_number()) %>% 
  group_by(multipolygon_id, polygon_id, linestring_id) %>% 
   arrange(row) %>% 
  summarize(area = polygon_area(cbind(x, y))) %>% 
  ungroup() #%>% 
  area %>%  mutate(area = area, 
                   hole = c(-1, 1)[(pmax(c(1, diff(multipolygon_id)), c(1, diff(polygon_id))) == 1) + 1])
}
xsf <- silicate::inlandwaters
df <- sfheaders::sf_to_df(xsf)
(mmarea <- calc_area(df))
#> `summarise()` has grouped output by 'multipolygon_id', 'polygon_id'. You can
#> override using the `.groups` argument.
#> # A tibble: 189 × 5
#>    multipolygon_id polygon_id linestring_id          area  hole
#>              <dbl>      <dbl>         <dbl>         <dbl> <dbl>
#>  1               1          1             1   2203380098.     1
#>  2               2          1             1    103039733.     1
#>  3               3          1             1 750145539384.     1
#>  4               3          1             2    100909770.    -1
#>  5               3          1             3   2203380098.    -1
#>  6               3          1             4     45585186.    -1
#>  7               3          1             5     89780792.    -1
#>  8               3          1             6    110755146.    -1
#>  9               3          1             7     59287239.    -1
#> 10               3          1             8     66189369.    -1
#> # ℹ 179 more rows
sum(mmarea$area * mmarea$hole)
#> [1] 1.923706e+12
sum(sf::st_area(xsf))
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
#> 1.923706e+12 [m^2]

Straightforward?

Well - ha ha - we will need some ease of use here, and some speed ups - but we aren’t restricted to one in-memory format, we have the area of every component part to use as we wish (the grouping tells us about holes and multipolygons), we can tell the orientation of the polygon or triangle and get its area, and when we have the right tools we’ll only be using R and C++. Oh and we can calculate with triangles stored compactly.

(a <- triangle_area(mm_tri$P[t(mm_tri$T), ]))
#>  [1] 0.03000 0.06900 0.03000 0.02000 0.05000 0.10000 0.05500 0.03000 0.03000
#> [10] 0.03000 0.03000 0.14550 0.05145 0.04500 0.11250 0.06250 0.06800
sum(a)
#> [1] 0.95895

## that is the area of the silicate::minimal_mesh (but the hole is included)

The key motivation here is flexibility and working with triangles.

Performance is compelling, compare silicate::tri_area() to this C++ version

library(silicate)
#> 
#> Attaching package: 'silicate'
#> The following objects are masked _by_ '.GlobalEnv':
#> 
#>     inlandwaters, minimal_mesh
#> The following object is masked from 'package:stats':
#> 
#>     filter
tas_tri <- RTriangle::triangulate(tas)
trxy <- tas_tri$P[t(tas_tri$T), ]
dim(trxy)
#> [1] 18291     2
rbenchmark::benchmark(R = {a <- tri_area(trxy)}, cpp = {b <- triangle_area(trxy)})
#>   test replications elapsed relative user.self sys.self user.child sys.child
#> 2  cpp          100   0.047    1.000     0.035    0.012          0         0
#> 1    R          100   0.215    4.574     0.202    0.012          0         0
sum(abs(a - b))
#> [1] 0.0004653583
mean(c(sum(a), sum(b)))
#> [1] 61672735794

---

Code of Conduct

Please note that the area project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.