voronoi function drops attributes from polygons

[Alternikaner]

Extending the example from built-in help to add arbitrary "readings" applicable to the polys:

library(terra)
#> terra version 1.2.10
wkt <- c("MULTIPOLYGON ( ((40 40, 20 45, 45 30, 40 40)), 
  ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35),(30 20, 20 15, 20 25, 30 20)))", 
         "POLYGON ((0 -5, 10 0, 10 -10, 0 -5))")
x <- vect(wkt)
x$readings <- rnorm(2)
x
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 2, 1  (geometries, attributes)
#>  extent      : 0, 45, -10, 45  (xmin, xmax, ymin, ymax)
#>  coord. ref. :  
#>  names       : readings
#>  type        :    <num>
#>  values      :   -1.524
#>                 -0.4625

v <- voronoi(x)
v
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 14, 0  (geometries, attributes)
#>  extent      : -55, 100, -65, 100  (xmin, xmax, ymin, ymax)
#>  coord. ref. :

^{Created on 2021-06-01 by the reprex package (v2.0.0)}

Is this intended behaviour? (Due to the increase in geometries, perhaps?) If so, how might one deal with polygon-based data for interpolation?

[rhijmans]

To keep the attributes, you can do:

v <- voronoi(as.points(x))

I suppose it would be reasonable to do this automatically for lines or polygons.

[Alternikaner]

Hi @rhijmans

Thanks for the lightning-quick response.

My 10 000-odd polys get converted to 1.3 million points, which causes the function to hang - or at least appear to. I've tried "st_simplify"-ing the original sf object before converting to a SpatVector, but this doesn't seem to have any effect on the number of points.

The voronoi function from the dismo package (which i believe you're also involved with -- many thanks for your numerous contributions!) doesn't have any issues dealing with my data, so if I stick with the raster+dismo combo, I do still have a working solution. Just wasn't sure if I was doing something wrong or if t-Voronoi treats attribute-polys differently than d-Voronoi by design.

[Nowosad]

Would not it be suitable here to use st_centroid() before calculating voronoi polygons?

[Alternikaner]

My data were collected from irregular adjoining areas of wildly different sizes, and is areal rather than point-based in nature.

Assuming voronoi generation takes the shapes of input polygons into account (i.e. doesn't simply calc centroids and work with a single point per poly in any case!) it is definitely the way I would wish to go in this particular case.

[rhijmans]

I have added the dismo/deldir method as an alternative algorithm. Now you can do:

library(terra)
f <- system.file("ex/lux.shp", package="terra")
x <- vect(f)

system.time(v1 <- voronoi(x, deldir=FALSE))
#   user  system elapsed 
#   0.16    0.02    0.18 
system.time(v2 <- voronoi(x, deldir=TRUE))
#   user  system elapsed 
#   1.16    0.00    1.16 

# with dismo
y <- as(as.points(x), "Spatial")
system.time(v3 <- dismo::voronoi(y))
#   user  system elapsed 
#   1.20    0.01    1.22 

So using deldir is much slower (because of the need for an lapply in R), but if it does not choke on large datasets, it is nice to have it (and I want to wean you off raster/dismo).

If you wanted the Voronoi diagram of the polygon centroids, you could do

z <- centroids(x)
v <- voronoi(z)

[Alternikaner]

Consider me weaned! :-)

Thanks so much for taking the time to add this, Robert.

As an aside, I was leaning on the rspatial.org interpolation tutorial for my switch-over attempt. It may be useful to point out the poly option here once the changes are on CRAN.

[Rohit-Satyam]

I am facing a similar issue. The above suggested method works if I do not use bnd argument but I have to use map boundary.

library(rnaturalearth)
library(mapview)
## Getting map
map <- ne_countries(country = "Australia", returnclass = "sf")

## mosquitoes 

mosquitoes<-read.csv("mosquitoes.csv")
map2 <- terra::vect(map)
x <- terra::vect(mosquitoes, geom=c("x","y"))
terra::crs(x) <- terra::crs(map2)
v <- voronoi(x,bnd=map2)
v

 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 1515, 0  (geometries, attributes)
 extent      : 113.339, 153.5695, -43.6346, -10.66819  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 

mosquitoes.zip