CODEWITHSUNDEEP
rbufferterra
Tarmo Remmel
Tarmo Remmel

Reputation: 11

terra::buffer() questionable functionality

When I try to create a circular polygon (using terra::buffer) to obtain a polygon with a specified area, the result is far from expected.

The code below illustrates the problem. Running R 4.2.2 on Windows (x86_64=w64-mingw32), using terra_1.7-3

Thanks.

I want to produce a circular polygon with a specified area. This is conceptually simple since the radius (r) is equal to the square-root of the Area (A) over Pi. r=sqrt(A/pi)

I then use the terra::buffer with DIST set as the computed radius to create the circular polygon, centered on a centroid point of a polygon from which I initially took the desired area. When I compute the terra::expanse() of this new circular polygon, it is nowhere near the desired area. Is this a bug or am I completely missing something? The desired area is ~1.7 larger than for the buffered circle.

library(terra)

# CONTROL WHICH POLYGON TO PROCESS WITHIN SpatVect
whichpoly <- 8

# SET GRAPHICAL ENVIRONMENT
plot.new()
par(mfrow=c(1,2))

# IMPORT DEMO SHAPEFILE AND PROJECT PROPERLY
v <- vect(system.file("ex/lux.shp", package="terra"))
test <- project(v, "EPSG:2150")

# PRODUCE STARTING MAP
plot(test, col="LightGrey", alpha=0.3)
plot(test[whichpoly], col="LightGreen", alpha=0.7, add=TRUE)
plot(centroids(test), add=TRUE)
plot(centroids(test[whichpoly]), col="Red", add=TRUE)
title(paste("Full Map: Centroids & Circle [", whichpoly, "]", sep=""))

# RATHER THAN BUILD A CIRCLE WITH THE EXTENT OF THE SHAPE, MATCH THE AREA TO THE SHAPE - NOT WORKING
# RADIUS OF A CIRCLE WITH SAME AREA AS SHAPE
# RADIUS = SQUARE-ROOT(AREA/PI)
radius2 <- sqrt(expanse(test[whichpoly])/pi)
buff2 <- buffer(centroids(test[whichpoly]), radius2)
plot(test[whichpoly])
plot(buff2, add=TRUE, col="cornflowerblue", alpha=0.5)

# AREAS
terra::expanse(test[whichpoly])
terra::expanse(buff2)
terra::expanse(test[whichpoly]) / terra::expanse(buff2)

Upvotes: 1

Views: 386

Answers (2)

margusl
margusl

Reputation: 17544

I'm afraid the used CRS (EPSG:2150) does not play well with that specific Shapefile (lux.shp, Luxembourg), it doesn't just get areas wrong (AREA attribute for 8th polygon, Grevenmacher, is 210 km^2) but also switches easting / northing and ends up having values out of UTM range. For projected CRS we could use +proj=utm +zone=31 as in the example bellow, though ?terra::expanse actually advises not to use planar coordinates for vectors:

For vector data, the best way to compute area is to use the longitude/latitude CRS. This is contrary to (erroneous) popular belief that suggest that you should use a planar coordinate reference system. This is done automatically, if transform=TRUE.

Long/lat SpatVector would also fine for terra::buffer()

library(terra)
#> terra 1.7.39

library(sf)
library(mapview)

whichpoly <- 8

v <- vect(system.file("ex/lux.shp", package="terra"))
test <- project(v, "+proj=utm +zone=31")

s <- st_as_sf(v)
mapviewOptions(basemaps = "CartoDB.Positron")
mapview(s, col.regions = "LightGreen") + 
  mapview(st_centroid(s)) + 
  mapview(st_centroid(s)[whichpoly,], col.regions = "red", alpha = 1)



radius2 <- sqrt(expanse(test[whichpoly])/pi)
buff2 <- buffer(centroids(test[whichpoly]), radius2)
plot(test[whichpoly])
plot(buff2, add=TRUE, col="cornflowerblue", alpha=0.5)



terra::expanse(test[whichpoly])
#> [1] 210354494
terra::expanse(buff2)
#> [1] 209358272
terra::expanse(test[whichpoly]) / terra::expanse(buff2)
#> [1] 1.004758

Created on 2023-08-03 with reprex v2.0.2

Upvotes: 1

Allan Cameron
Allan Cameron

Reputation: 174468

The discrepancy vanishes if you use transform = FALSE inside extent

library(terra)
#> terra 1.7.3

whichpoly <- 8

plot.new()
par(mfrow=c(1,2))

v <- vect(system.file("ex/lux.shp", package="terra"))
test <- project(v, "EPSG:2150")

plot(test, col="LightGrey", alpha=0.3)
plot(test[whichpoly], col="LightGreen", alpha=0.7, add=TRUE)
plot(centroids(test), add=TRUE)
plot(centroids(test[whichpoly]), col="Red", add=TRUE)
title(paste("Full Map: Centroids & Circle [", whichpoly, "]", sep=""))

radius2 <- sqrt(expanse(test[whichpoly], transform = FALSE)/pi)
buff2 <- buffer(centroids(test[whichpoly]), radius2)
plot(test[whichpoly])
plot(buff2, add=TRUE, col="cornflowerblue", alpha=0.5)



# AREAS
terra::expanse(test[whichpoly], transform = FALSE)
#> [1] 360211556
terra::expanse(buff2, transform = FALSE)
#> [1] 358732072
terra::expanse(test[whichpoly]) / terra::expanse(buff2)
#> [1] 1.004125

Created on 2023-08-02 with reprex v2.0.2

Upvotes: 1

Related Questions