I'm starting from the attached raster file which is a world map of 5 climate zones. I then convert it into a simple feature. My goal is to get rid of those small polygons contained in the larger ones so that I'm only visualizing the outlines of the larger patches of the five different zones. As a bonus it would be nice to make those line smoother because the edges are not visually pleasing but being a categorical quantity I can't really interpolate the values.
library(stars)
x=read_stars('climi.tif')
#convert into sf
tes=st_as_sf(x,as_points=FALSE,merge=TRUE,use_integer=TRUE)
#load a simple background map
map.background=sf::st_as_sf(rworldmap::countriesCoarseLessIslands)
ggplot()+ geom_sf(data=map.background,size=0.1) +geom_sf(data=tes,aes(geometry=geometry,color=factor(climi.tif)),fill=NA,size=0.7)
I get this, but as you can see there are many small polygons of a few pixels inside very large polygons, this makes the visualization confusing. I would like to merge those small polygons with the larger ones that surround them.

I have tried to calculate area and remove the smaller polygons and then use sf_remove_holes, which improves the result but I still get too many smaller polygons or NA areas which should be colored.
tes$area= st_area(tes);tes$area=as.numeric(tes$area)
subtes=subset(tes,tes$area>199426036400)
subtes=sfheaders::sf_remove_holes(subtes)
But the some small polygons are still visible and additionally it seems that some areas don't seem to belong to any zone. I have tried also st_contains and raster::clump() but without success. Clump for example removes one of the five categories. I would like to form continuous areas that are easy to distinguish between each other, even if I have to draw double borders.
