4

I created a spatial polygon file from a classified gridded dataset using the grid2polygons function in R. I would now like to take that polygon file and merge any polygons containing less than 6 pixels with its nearest largest neighbor. In the image link below, you can see an example of small polygons circled in blue. I would like those merged with the larger surround polygon. Alternatively, if there is a way to create a minimum threshold of 6 clumped pixels prior to the conversion to polygon stage, any suggestions for that method would be great as well.

Polygon Image

There is a function in Gdal: gdal_sieve.py that can do exactly this (http://www.gdal.org/gdal_sieve.html) if you would like to get an idea of exactly what I would like to do. Does anyone know of a way to do this in R?

2

You can call gdal_sieve using system2() in R if you have your system environment variables set up correctly, e.g.:

system2('gdal_sieve', args = c('-st', 6, infile.tif, outfile.tif))

then vectorise the output with

system2('gdal_polygonize', args = c(outfile.tif, '-f', 'GPKG', 'outfile.gpkg'))

You can do all this using GRASS via R too. The process is a bit more involved, but you have more control over some aspects, and it handles no-data areas properly. I've blogged it here but to summarise:

  • read your raster into GRASS
  • use r.to.vect to vectorise the raster
  • use v.clean with tool 'rmarea' and an appropriate threshold to remove extraneous tiny polygons
  • if you want something with smoother boundaries, convert the results of the above process back into a raster, then convert it to vector again using the flag 's'.

Be aware that the topology in the resulting outputs can be problematic and that no easy fixes appear to exist (gosh, I'd be happy to be wrong about this). So the outputs can be fine for interactive use, but subsequent geometric operations like intersect and dissolve can fail unless you make some manual edits.

  • Thank you for this information! I can use this for some of my other work. However, the problem I have is that I have been tasked to complete this entirely in R which is why I was hoping there would be some code that might do this. As I have yet to find an answer to this anywhere on the interweb, this may be one of those times I get motivated to learn how to edit function source code. – Rachel Jun 6 '17 at 3:37
  • Ah, ok. R spatial packages rely heavily on external libraries for geometric ops anyway (GEOS mostly, and some JS libs like mapshaper) so that might not have been a reasonable task to set... – obrl_soil Jun 6 '17 at 9:51
  • Yeah. I realized this might be an issue and expressed that to my boss. I figured if there wasn't an option in either rgdal or gdalUtils, it wouldn't exist without installing and calling gdal or writing functions for a new package. – Rachel Jun 7 '17 at 0:53

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.