drop-large-triangles.R

library(sf)

## Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.3

east <- c(438858.0, 438865.0, 438872.0, 438879.1, 438886.1, 438893.1, 438900.2, 438907.3, 438914.4,
          438921.3, 438928.2, 438935.1, 438941.8, 438942.5, 438936.2, 438929.8, 438923.4, 438916.9,
          438910.3, 438903.7, 438897.2, 438890.9, 438884.9, 438879.3, 438873.3, 438867.2, 438860.7,
          438864.1, 438870.9, 438877.6, 438931.4, 438938.4, 438945.3, 438951.9, 438945.4, 438938.6,
          438931.8, 438878.9, 438871.8, 438865.0, 438869.1, 438875.6, 438882.2, 438936.0, 438942.6,
          438948.4, 438953.2, 438954.4, 438948.3, 438942.0, 438935.6, 438929.0, 438922.4, 438915.7,
          438909.0, 438902.6, 438896.4, 438890.8, 438885.5, 438880.0, 438874.0, 438867.9, 438869.8,
          438876.4, 438882.9, 438889.4, 438895.8, 438902.2, 438908.6, 438914.9, 438920.9, 438927.0,
          438933.1, 438939.4, 438945.9, 438952.5)
north <- c(5991835, 5991834, 5991834, 5991833, 5991832, 5991831, 5991830, 5991829, 5991829, 5991828,
           5991827, 5991826, 5991825, 5991843, 5991843, 5991844, 5991845, 5991846, 5991847, 5991848,
           5991848, 5991849, 5991850, 5991850, 5991851, 5991852, 5991853, 5991873, 5991872, 5991871,
           5991864, 5991863, 5991862, 5991879, 5991880, 5991881, 5991881, 5991888, 5991889, 5991890,
           5991906, 5991905, 5991905, 5991898, 5991897, 5991896, 5991895, 5991916, 5991917, 5991918,
           5991919, 5991920, 5991920, 5991921, 5991922, 5991923, 5991924, 5991924, 5991925, 5991926,
           5991926, 5991927, 5991945, 5991944, 5991944, 5991943, 5991942, 5991941, 5991940, 5991940,
           5991939, 5991938, 5991937, 5991936, 5991936, 5991935)
df <- data.frame(x = east, y = north)
p1 <- st_as_sf(x = df, coords = c("x", "y"), remove = F,
               crs = "+proj=utm +zone=20 +south +datum=WGS84 +units=m +no_defs")

## an alternative approach to this question: 
## https://gis.stackexchange.com/questions/254836/sf-lines-to-polygons-with-holes
library(sfdct)
p <- st_cast(ct_triangulate(p1))

## all POINT, returning one feature triangulated

## Warning in st_cast.sf(ct_triangulate(p1)): repeating attributes for all
## sub-geometries for which they may not be constant

a <- as.numeric(st_area(p$geometry))
## use a finite element approach and drop large triangles
pp <- st_union(p[a < quantile(a, 0.95), ])
plot(pp, col = "grey")
plot(p1, add = TRUE)