Fast way to convert raster to polygon shapefile in R

Question

I have a single-band raster showing different land cover classes. I want to convert it to a polygon shapefile. I tried to use the rasterToPolygons function with the option to dissolve neighboring similar values. It takes hours to complete it. Are there any alternative fast options? Here is additonal information about the file

class      : RasterLayer 
:Formal class '.RasterFile' [package "raster"] with 13 slots
dimensions : 618, 1543, 953574  (nrow, ncol, ncell)
resolution : 0.001619785, 0.001619808  (x, y)
extent     : -107.9999, -105.5005, 48.99976, 50.00081  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : testNALCMSlatlon252 
values     : 1, 18  (min, max)

Answer 1

Use the terra package instead of raster and use as.polygons to polygonize.

Sample data maker:

library(terra)
maker = function(x,y){rast(matrix(sample(1:18, x*y, TRUE), x, y))}

# same scale as yours:
r18 = maker(600,1500)

this creates a random raster which is probably a worse-case scenario for polygonising - lots of tiny single-cell polygons.

How long to convert to polygons...

> system.time({p18 = as.polygons(r18)})
   user  system elapsed                   
  3.866   0.024   3.890 

okay 4 seconds. Quick enough? I suspect it works fast because it knows about grid adjacency, whereas raster's use of rasterToPolygons uses rgeos which has to compute adjacency from the polygon coordinates which is sloooow...

So what have we got?

> p18
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 18, 1  (geometries, attributes)
 extent      : 0, 1500, 0, 600  (xmin, xmax, ymin, ymax)
 coord. ref. :  
 names       : lyr.1
 type        : <int>
 values      :     1
                   2
                   3

Its a SpatVector object with 18 features, so one feature for all polygons in each class, as a multi-polygon feature.

Let's illustrate with a small example with fewer classes:

> r2 = round(maker(10,10)/6)
> plot(r2)                                

> p2 = as.polygons(r2)
> plot(p2)                                
> plot(p2, col=1:4)

And you can see the polygons matching the raster categories.

If you're not familiar with SpatVector classes then neither am I, but I'm sure you can do stuff with them like convert to sf classes to use functions from the sf package etc.