0

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.

3
  • The link for the data on Google Drive is locked. Please open the permissions or remove the link. Nov 5 '21 at 20:45
  • As per the help center please do not include chit chat like thanks and signatures in your posts.
    – PolyGeo
    Nov 6 '21 at 21:14
  • I have updated the link to the data, to make it full permissions, as far as I can tell. Thanks. Nov 7 '21 at 19:26
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')

3
  • Thanks, I will try this. I also updated data link above. Nov 7 '21 at 19:29
  • I've tried this now and had two issues: 1) I get an error when running st_as_sf() on a SpatVector ("Error in UseMethod("st_as_sf") : no applicable method for 'st_as_sf' applied to an object of class "SpatVector"). 2) Even just using the as.polygons function, I still get the three meter offset. So, for the latter, I must be not understanding something about how extents are maintained during the (generic) raster to vector conversion. Nov 8 '21 at 13:28
  • Update: confirmed this offset also happens with normal raster functions (raster::rasterToPolygons). So then the question is, how does one avoid this shift? For high-resolution data it is problematic for sure. Nov 8 '21 at 13:44
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.

6
  • These are zonal polygons in a small landscape, there are only a couple thousand and it's nice to have them as vector outlines to overlay on other data. That said, what is the preferred "Zonal Statistics as Table" workflow in R? It seemed like using raster::extract or exactextractR were the most straightforward ways to get raster data into a data frame. Nov 16 '21 at 18:59
  • The raster/terra approach would be to use zonal. That avoids the need to create polygons, but you must have matching rasters. For extracting I would use terra::extract over raster::extract; or else exactextracr Nov 16 '21 at 19:05
  • Running your code worked fine. However, even when I start with your raster as template and insert only my own pixel values, the issue comes back. I believe now that the issue is with NA values. When I set the NA values to 0, the issue goes away. Not sure if this is specific to my dataset or a general issue. Thanks! Nov 16 '21 at 19:39
  • If you want to keep the NAs you can do as.polygons(x, na.rm=FALSE) Nov 16 '21 at 21:00
  • Great. I think this is the actual answer to my original question. But at the end of the day, this question, as posed, is probably not that useful to others. A bit of a red herring. Nov 17 '21 at 13:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.