rgdalterra
Inner buffer in terra giving higher area than its template
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
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
#;-)
#;-) ──────────────────────────────────────────────────────────────────────────────
Solution
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.
- Floating point less-than-equal comparisons after addition and subtraction
- ggplot fill spilling over shape when used with different geoms
- Compute mutual gaze from individual gazes
- Inverse of gamma function with R
- Conditionally mutating 2 dataframes
- Easiest way to save an S4 class
- Updating values of three dot ellipsis in R
- How to make a multifactor model in pROC?
- Using hidden API to scrape website with JavaScript pager
- Map colours not showing correctly, all SpatRaster values converted to 1 when using crop and mask
- Hierarchical Bayes for R or Python
- How to insert the "Copy to clipboard" icon in an r markdown document
- Name new variables with purrr::map()
- Show every time stamp in plotly x-axis
- Represent the selectInput value in an infoBox
- Display the selectInput value in a R shiny widget box
- SelectInput Value update based on previous SelectInput in R shiny
- Using Conditional Panel with PickerInput in R
- Display multiple infoboxes within one reactive function
- R shiny sever problems - reactive global environment
- Show / Hide selectinput in R shiny based on conditions
- Update data from selectInput
- Usage of UIOutput in multiple menuItems
- Automating the slider of SelectInput
- ggplot - geom_point: setting manual fill color for point jitter without changing boxplot fill color
- why does assigning to the ellipsis `... <- 1` not throw an error?
- checkboxGroupInput delete table output in a shiny app
- ShinyApp: How to pass plot from one module into another one in order to download it in a report?
- Names of a numeric vector with HTML
- How to place elements above the dashboardHeader