st_as_sf() of stars object reduces the object extent

Question

I am attempting to do a raster to polygon conversion using the stars package. The workflow i've seen elsewhere (and have used successfully previously) is on my "zones" raster is:

zones_stars <- st_as_stars(zones) %>% st_as_sf(merge = TRUE)

When I break this pipeline apart for diagnostics, the first step (st_as_stars) seems to result in maintenance of original object extents in raster format.

(Original) Zones raster:

class      : RasterLayer 
dimensions : 2744, 5907, 16208808  (nrow, ncol, ncell)
resolution : 0.02, 0.02  (x, y)
extent     : 488532.1, 488650.3, 4424048, 4424103  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=NAD83 +units=m +no_defs 
source     : memory
names      : layer 
values     : 172, 12263  (min, max)

Output of st_as_stars:

attribute(s), summary of first 1e+05 cells:
       Min. 1st Qu. Median Mean 3rd Qu. Max.  NA's
layer    NA      NA     NA  NaN      NA   NA 1e+05
dimension(s):
  from   to  offset delta                       refsys point values x/y
x    1 5907  488532  0.02 +proj=utm +zone=18 +datum...    NA   NULL [x]
y    1 2744 4424103 -0.02 +proj=utm +zone=18 +datum...    NA   NULL [y]

So far so good. The second step (st_as_sf), converts to vector polygons but introduces a several meter offset:

geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 488532.1 ymin: 4424051 xmax: 488650.3 ymax: 4424102
CRS:            +proj=utm +zone=18 +datum=NAD83 +units=m +no_defs
First 10 features:
   layer                       geometry
1    332 POLYGON ((488648.2 4424102,...
2    335 POLYGON ((488647.7 4424102,...
3    172 POLYGON ((488646.1 4424101,...
4    338 POLYGON ((488646.9 4424101,...

Note particularly, the ymin value shift from original to final output. Any thoughts on this?

Here is the original zones data, if you'd like to try to re-create: https://bit.ly/31yPfhz

EDIT: The issue appears to be related to the inclusion of NA values in the raster. If you convert those to a number (e.g., 0) then there is no issue. I do not know if this applies more generally or is data-specific.

Answer 1

If you create polygons form a raster with values, the NA's are removed by terra::as.polygons unless you set na.rm=FALSE. If the NAs are at the edges (as they often are), the extent can get reduced.

Otherwise, I do not see a change in the extent:

library(terra)
r <- rast(xmin=488532.1, xmax=488650.3, ymin=4424048, ymax=4424103, nrow=2744, ncol=5907)
x <- as.polygons(r)

r
#class       : SpatRaster 
#dimensions  : 2744, 5907, 1  (nrow, ncol, nlyr)
#resolution  : 0.02001016, 0.02004373  (x, y)
#extent      : 488532.1, 488650.3, 4424048, 4424103  (xmin, xmax, ymin, ymax)
coord. ref. :  

x
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 16208808, 0  (geometries, attributes)
# extent      : 488532.1, 488650.3, 4424048, 4424103  (xmin, xmax, ymin, ymax)
# coord. ref. :  

And then to sf, the extent (bbox) is still the same (only the order of appearance is different)

s <- sf::st_as_sf(x)
sf::st_bbox(s)
# xmin      ymin      xmax      ymax 
# 488532.1 4424048.0  488650.3 4424103.0 

The message no applicable method for 'st_as_sf' applied to an object of class "SpatVector" suggests that you are using an antique version of sf, and that you should update your packages.

I do wonder why you would create polygons from a large raster --- that is almost always unnecessary and extremely inefficient.

Answer 2

If you want to try a different approach I would recommend using terra (Rcpp based replacement for raster package). To read a raster it would be r <- terra::rast("x") and, to convert to polygon you use terra:as.polygons.

Since this results in a Vect object class (which are a pain to work with outside of terra) you can nest the polygon conversion in a coercion to sf

p <- sf::st_as_sf(as.polygons(x, trunc=TRUE, dissolve=TRUE))

This will be a MULTIPOLYGON geometry so, to explode to single part, you can cast the polygons p <- sf::st_cast(p, 'POLYGON')