Inner buffer in terra giving higher area than its template

Question

Strange issue with inner buffer in terra

Goal

My goal is to make 2 shapefiles 1 with the strip of land 100 kms from the coast and another one with all the world but that. In order to do that I am using the following packages:

library(terra)
#;-) terra 1.6.17
library(rworldxtra)
#;-) Loading required package: sp
library(tidyterra)
#;-) 
#;-) Attaching package: 'tidyterra'
#;-) The following object is masked from 'package:stats':
#;-) 
#;-)     filter
library(tidyverse)

I first use the countriesHigh data base from the rworldxtra package

Then I transform it to a spatVector with terra and create a new variable to agreggate over with

World <- terra::vect(countriesHigh)

World$World <- "Yes"

Then I select only that variable, aggregate by it and make it valid

World <- World[, "World"]

World <- terra::makeValid(World)

World <- terra::aggregate(World, by = "World", cores = 3)

I fill the holes to fix problems and make it valid again, and since aggregate generates a new layer, I reselect only the one I intended to use

World <- terra::fillHoles(World) %>% terra::project("+proj=moll")

World <- terra::makeValid(World)

World <- World[, "World"]

I changed the projection to mollweide to avoid dateline problems This leads me to the following map:

ggplot() +
  geom_spatvector(data = World, fill = "green") +
  theme_bw() +
  ggtitle("World")

I also calculate the area in SqKms

WorldArea <- terra::expanse(World, unit = "km") %>% sum()
#;-) Warning in x@ptr$area(unit, transform, double()): GDAL Error 1: Point outside of
#;-) projection domain

#;-) Warning in x@ptr$area(unit, transform, double()): GDAL Error 1: Point outside of
#;-) projection domain

#;-) Warning in x@ptr$area(unit, transform, double()): GDAL Error 1: Point outside of
#;-) projection domain

#;-) Warning in x@ptr$area(unit, transform, double()): GDAL Error 1: Point outside of
#;-) projection domain

The World area is 67,398,716 sq km

Then I generate an inner world, that is one with a negative buffer of 100 kms

World_Inner <- terra::buffer(World, width = -100000)

World_Inner <- terra::makeValid(World_Inner)

As you can see this looks fine

ggplot() +
  geom_spatvector(data = World_Inner, fill = "green") +
  theme_bw() +
  ggtitle("Inner")

Areas close to the coast like Panama and Denmark disappear, however when I check the area:

InnerArea <- terra::expanse(World_Inner, unit = "km") %>% sum()

The Inner area is 114,840,316, that is r (terra::expanse(World_Inner, unit = “km”)/terra::expanse(World, unit = “km”))*100, 2) percent of the world so bigger than the world?

That is particuallarly weird when I see that it actually fits inside and looks fine

ggplot() +
  geom_spatvector(data = World, fill = "green") +
  geom_spatvector(data = World_Inner, fill = "red") +
  theme_bw() +
  ggtitle("Inner")

Finally I generate the coastal area

Coastal_World <- erase(World, World_Inner)

Which looks like this:

ggplot() +
  geom_spatvector(data = Coastal_World, fill = "green") +
  theme_bw() +
  ggtitle("Coastal")

again fine, although the area is bigger than I expected

CoastalArea <- terra::expanse(Coastal_World, unit = "km")
#;-) Warning in x@ptr$area(unit, transform, double()): GDAL Error 1: Point outside of
#;-) projection domain

#;-) Warning in x@ptr$area(unit, transform, double()): GDAL Error 1: Point outside of
#;-) projection domain

#;-) Warning in x@ptr$area(unit, transform, double()): GDAL Error 1: Point outside of
#;-) projection domain

#;-) Warning in x@ptr$area(unit, transform, double()): GDAL Error 1: Point outside of
#;-) projection domain

The Costal area is 21,134,419 sq km, that is 31.36 percent of the world.

Any idea what I am doing wrong? The only issue I can see are this warnings about points being outside the projection domain, I cant figure out that issue

The generated shapefiles are here in case someone wants to check them

^{Created on 2022-11-14 by the reprex package (v2.0.1)}

Standard output and standard error

ℹ The Stack Overflow venue "so" is an alias for the default GitHub venue "gh".
There is no need to specify the venue.
Error in (function (x)  : attempt to apply non-function

Session info

sessioninfo::session_info()
#;-) ─ Session info ───────────────────────────────────────────────────────────────
#;-)  setting  value
#;-)  version  R version 4.2.2 Patched (2022-11-10 r83330)
#;-)  os       Ubuntu 20.04.5 LTS
#;-)  system   x86_64, linux-gnu
#;-)  ui       X11
#;-)  language en_US:en
#;-)  collate  en_US.UTF-8
#;-)  ctype    en_US.UTF-8
#;-)  tz       Europe/Copenhagen
#;-)  date     2022-11-14
#;-)  pandoc   2.17.1.1 @ /usr/lib/rstudio/bin/quarto/bin/ (via rmarkdown)
#;-) 
#;-) ─ Packages ───────────────────────────────────────────────────────────────────
#;-)  package     * version date (UTC) lib source
#;-)  assertthat    0.2.1   2019-03-21 [1] CRAN (R 4.2.0)
#;-)  backports     1.4.1   2021-12-13 [1] CRAN (R 4.2.0)
#;-)  broom         1.0.1   2022-08-29 [1] CRAN (R 4.2.1)
#;-)  cellranger    1.1.0   2016-07-27 [1] CRAN (R 4.2.0)
#;-)  class         7.3-20  2022-01-13 [4] CRAN (R 4.1.2)
#;-)  classInt      0.4-8   2022-09-29 [1] CRAN (R 4.2.2)
#;-)  cli           3.4.1   2022-09-23 [1] CRAN (R 4.2.1)
#;-)  codetools     0.2-18  2020-11-04 [4] CRAN (R 4.0.3)
#;-)  colorspace    2.0-3   2022-02-21 [3] CRAN (R 4.1.2)
#;-)  crayon        1.5.2   2022-09-29 [3] CRAN (R 4.2.1)
#;-)  curl          4.3.3   2022-10-06 [3] CRAN (R 4.2.1)
#;-)  DBI           1.1.3   2022-06-18 [1] CRAN (R 4.2.1)
#;-)  dbplyr        2.2.0   2022-06-05 [1] CRAN (R 4.2.0)
#;-)  digest        0.6.30  2022-10-18 [1] CRAN (R 4.2.1)
#;-)  dplyr       * 1.0.10  2022-09-01 [1] CRAN (R 4.2.1)
#;-)  e1071         1.7-12  2022-10-24 [1] CRAN (R 4.2.2)
#;-)  ellipsis      0.3.2   2021-04-29 [3] CRAN (R 4.0.5)
#;-)  evaluate      0.18    2022-11-07 [3] CRAN (R 4.2.2)
#;-)  fansi         1.0.3   2022-03-24 [3] CRAN (R 4.1.3)
#;-)  farver        2.1.1   2022-07-06 [3] CRAN (R 4.2.1)
#;-)  fastmap       1.1.0   2021-01-25 [3] CRAN (R 4.0.3)
#;-)  forcats     * 0.5.2   2022-08-19 [1] CRAN (R 4.2.1)
#;-)  fs            1.5.2   2021-12-08 [3] CRAN (R 4.1.2)
#;-)  generics      0.1.3   2022-07-05 [1] CRAN (R 4.2.1)
#;-)  ggplot2     * 3.4.0   2022-11-04 [1] CRAN (R 4.2.2)
#;-)  glue          1.6.2   2022-02-24 [3] CRAN (R 4.1.2)
#;-)  gtable        0.3.1   2022-09-01 [3] CRAN (R 4.2.1)
#;-)  haven         2.5.0   2022-04-15 [1] CRAN (R 4.2.0)
#;-)  highr         0.9     2021-04-16 [3] CRAN (R 4.0.5)
#;-)  hms           1.1.2   2022-08-19 [1] CRAN (R 4.2.1)
#;-)  htmltools     0.5.3   2022-07-18 [3] CRAN (R 4.2.1)
#;-)  httr          1.4.4   2022-08-17 [1] CRAN (R 4.2.1)
#;-)  jsonlite      1.8.3   2022-10-21 [3] CRAN (R 4.2.1)
#;-)  KernSmooth    2.23-20 2021-05-03 [4] CRAN (R 4.0.5)
#;-)  knitr         1.40    2022-08-24 [1] CRAN (R 4.2.1)
#;-)  lattice       0.20-45 2021-09-22 [4] CRAN (R 4.2.0)
#;-)  lifecycle     1.0.3   2022-10-07 [3] CRAN (R 4.2.1)
#;-)  lubridate     1.8.0   2021-10-07 [1] CRAN (R 4.2.0)
#;-)  magrittr      2.0.3   2022-03-30 [3] CRAN (R 4.1.3)
#;-)  mime          0.12    2021-09-28 [3] CRAN (R 4.1.1)
#;-)  modelr        0.1.8   2020-05-19 [1] CRAN (R 4.2.0)
#;-)  munsell       0.5.0   2018-06-12 [3] CRAN (R 4.0.0)
#;-)  pillar        1.8.1   2022-08-19 [1] CRAN (R 4.2.1)
#;-)  pkgconfig     2.0.3   2019-09-22 [3] CRAN (R 4.0.0)
#;-)  proxy         0.4-27  2022-06-09 [1] CRAN (R 4.2.1)
#;-)  purrr       * 0.3.5   2022-10-06 [3] CRAN (R 4.2.1)
#;-)  R.cache       0.16.0  2022-07-21 [1] CRAN (R 4.2.1)
#;-)  R.methodsS3   1.8.2   2022-06-13 [1] CRAN (R 4.2.1)
#;-)  R.oo          1.25.0  2022-06-12 [1] CRAN (R 4.2.1)
#;-)  R.utils       2.12.0  2022-06-28 [1] CRAN (R 4.2.1)
#;-)  R6            2.5.1   2021-08-19 [3] CRAN (R 4.1.1)
#;-)  Rcpp          1.0.9   2022-07-08 [3] CRAN (R 4.2.1)
#;-)  readr       * 2.1.3   2022-10-01 [1] CRAN (R 4.2.1)
#;-)  readxl        1.4.0   2022-03-28 [1] CRAN (R 4.2.0)
#;-)  rematch2      2.1.2   2020-05-01 [3] CRAN (R 4.0.0)
#;-)  reprex        2.0.1   2021-08-05 [1] CRAN (R 4.2.0)
#;-)  rlang         1.0.6   2022-09-24 [1] CRAN (R 4.2.1)
#;-)  rmarkdown     2.17    2022-10-07 [1] CRAN (R 4.2.1)
#;-)  rstudioapi    0.14    2022-08-22 [1] CRAN (R 4.2.1)
#;-)  rvest         1.0.2   2021-10-16 [1] CRAN (R 4.2.0)
#;-)  rworldxtra  * 1.01    2012-10-03 [1] CRAN (R 4.2.1)
#;-)  scales        1.2.1   2022-08-20 [3] CRAN (R 4.2.1)
#;-)  sessioninfo   1.2.2   2021-12-06 [3] CRAN (R 4.1.2)
#;-)  sf            1.0-9   2022-11-08 [1] CRAN (R 4.2.2)
#;-)  sp          * 1.5-1   2022-11-07 [1] CRAN (R 4.2.2)
#;-)  stringi       1.7.8   2022-07-11 [3] CRAN (R 4.2.1)
#;-)  stringr     * 1.4.1   2022-08-20 [1] CRAN (R 4.2.1)
#;-)  styler        1.7.0   2022-03-13 [1] CRAN (R 4.2.1)
#;-)  terra       * 1.6-17  2022-09-10 [1] CRAN (R 4.2.1)
#;-)  tibble      * 3.1.8   2022-07-22 [1] CRAN (R 4.2.1)
#;-)  tidyr       * 1.2.1   2022-09-08 [1] CRAN (R 4.2.1)
#;-)  tidyselect    1.2.0   2022-10-10 [1] CRAN (R 4.2.1)
#;-)  tidyterra   * 0.3.0   2022-10-12 [1] CRAN (R 4.2.1)
#;-)  tidyverse   * 1.3.1   2021-04-15 [1] CRAN (R 4.2.0)
#;-)  tzdb          0.3.0   2022-03-28 [1] CRAN (R 4.2.0)
#;-)  units         0.8-0   2022-02-05 [1] CRAN (R 4.2.0)
#;-)  utf8          1.2.2   2021-07-24 [3] CRAN (R 4.1.0)
#;-)  vctrs         0.5.0   2022-10-22 [1] CRAN (R 4.2.1)
#;-)  withr         2.5.0   2022-03-03 [3] CRAN (R 4.1.3)
#;-)  xfun          0.34    2022-10-18 [3] CRAN (R 4.2.1)
#;-)  xml2          1.3.3   2021-11-30 [3] CRAN (R 4.1.2)
#;-)  yaml          2.3.6   2022-10-18 [3] CRAN (R 4.2.1)
#;-) 
#;-)  [1] /home/au687614/R/x86_64-pc-linux-gnu-library/4.2
#;-)  [2] /usr/local/lib/R/site-library
#;-)  [3] /usr/lib/R/site-library
#;-)  [4] /usr/lib/R/library
#;-) 
#;-) ──────────────────────────────────────────────────────────────────────────────

https://stackoverflow.com/questions/68308376/why-do-terracellsize-and-rasterarea-produce-different-estimates-of-raste — Hansel Palencia, Nov 14 '22 at 13:58
I think it could be the fact that the terra calculates area differently? See link above... — Hansel Palencia, Nov 14 '22 at 14:13
Both areas are calculated by terra, maybe I am missing something @HanselPalencia, but I dont think that if both areas are calculated by terra there should be this kind of mistake' — Derek Corcoran, Nov 14 '22 at 14:14
I also just looked up the total area of the world which is 509 600 000 square km, while the total area of land in the world is 148 326 000 km. Quite a stark difference from your reported 68 million? — Hansel Palencia, Nov 14 '22 at 14:16
To be fair though if inner was 114,840,316 and world land is 148,326,000 which is around 77.4% — Hansel Palencia, Nov 14 '22 at 14:19
Or even better if we use the coastal value - 21,134,419 divided by the total land area value 148,326,000 that equals only 14.2% of all land area. Either way I'm looking at it, there seems to be a problem with your initial world area value. — Hansel Palencia, Nov 14 '22 at 14:20
Not addressing the problem, but best to periodically update via git clone to most recent development version and `terra`, now about `1.6.34` or so. And what version is you `GDAL`? — Chris, Nov 14 '22 at 14:50
Tried it with the new version, and still the same results @Chris — Derek Corcoran, Nov 14 '22 at 18:36

Robert Hijmans · Accepted Answer · 2022-11-15T03:27:36.733

The results are reasonable when doing this, that is, compute area ignoring distortion from the coordinate reference system.

library(terra)
library(geodata)
w <- world(path=".")
w <- project(w, "+proj=moll") |> aggregate()
aw <- expanse(w, unit = "km", transform=F)

b <- buffer(w, width = -100000)
ab <- expanse(b, unit = "km", transform=F)

ab / aw
#[1] 0.7890138

But things go wrong when using transform=TRUE. In that case, the polygons are projected to lon/lat before computing the area. That should be more precise, but something is going wrong.

aw2 <- expanse(w, unit = "km", transform=T)
ab2 <- expanse(b, unit = "km", transform=T)
ab2 / aw2
#[1] 2.11017

Because aw2 is much too small. This is related to the warnings you are getting. Illustrated here

d <- disagg(w)
afeuras <- d[24,]
expanse(afeuras)
#[1] 0
#Warning messages:
#1: Point outside of projection domain (GDAL error 1)

I need to investigate to see what can be done.

Thanks @Robert, I think I know why this is, when reprojecting to lon lat a lot of things go wrong specially with the dateline, so polygons get super distorted, that is why I transformed into moll in the first place, i tried with terra::normalize.longitude but it does not do the trick — Derek Corcoran, Nov 15 '22 at 15:10
PS Robert is any of the datasets in Geodata a good idea for doing this? — Derek Corcoran, Nov 15 '22 at 15:22

Inner buffer in terra giving higher area than its template

Strange issue with inner buffer in terra

Goal

1 Answers1