Search code examples
rk-meansspatialrasterr-sp

Automatically reclassify small SpatialPolygons inserted into large SpatialPolygons using R

I would like to assign small polygons nested in larger polygons the same values as those of larger polygons. In figure 1 you can see the small polygons in raster format:

Raster with small polygons within larger polygons

and in figure 2 in SpatialPolygons as individual polygons: Disaggregated polygons

These polygons are results of sorting by k-means, generating raster, and using the rasterFromXYZ function (code below):

mydata.26.raster <- rasterFromXYZ(as.data.frame(mydata.26.coord[,c("x", "y",       "cls_26.cluster")]),res=5, crs=crs)

and then rasterToPolygons function I was able to separate the polygons (code below):

zona.26.pol<- rasterToPolygons(zona.26.raster$cls_26.cluster,dissolve=TRUE)
zona.26.pol <- disaggregate(zona.26.pol)

Here's zona.26.pol if you want to help It is in .shp format.

And manually I reclassified the polygons and finally added them using the same classes. After manually assigning the values by me, the result that I would like to achieve automatically (creating rules) is in figure 3:

Reclassified and aggregated final image

Every help is welcome!

Solution

  • This will remove the small nested polygons based on their size alone and then remove the holes left in the larger, remaining polygons. This works for your example but may fail if you have larger nested polygons you are wanting to remove. In that case, we would have to figure out how to identify the geometries that are 'nested' 

    library(sf)
library(units)
library(nngeo)
min_polygon_area  <- 10000 #set minimum size of a nested polygon you would like to remove
units(min_polygon_area) <- as_units('m^2') #as defined below by st_area
zona.26.pol <- st_read(file.path(workDir, 'zona.26.pol.shp'))
st_crs(zona.26.pol) <-  '+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs' #define crs
zona.26.pol$area <- st_area(zona.26.pol)
zona.26.pol$area #note area is in m^2
plot(zona.26.pol[,'cls_26_']) #as-is plot

    enter image description here

    zona.26.pol <- zona.26.pol[zona.26.pol$area>min_polygon_area, ]
plot(zona.26.pol[,'cls_26_']) #small polygons removed; holes remaining

    enter image description here

    zona.26.pol_no_holes <- st_remove_holes(zona.26.pol)
plot(zona.26.pol_no_holes[,'cls_26_']) #holes removed

    enter image description here Note that I used the sf package to read-in the shapefile, in order to utilize the st_remove_holes function from the nngeo package, but I typically use raster and sp packages.