1

I have a large geotiff with NDVI differences (values -10'000 to +10'000), that I use to create polygons. My workflow currently is to polygonize a binary mask of it and then filter areas that are too small. The problem with this approach is, that there are a lot of tiny polygons that eventually I don't need, but that are created with gdal_polygonize (taking days of computation time).

Can I filter the binary mask before polygonizing, only creating polygons larger than AREA_THRESHOLD (e.g. 5 connected pixels = 500 m²) ?


My current workflow is:

  1. filter my raster for values > THRVALUE (I need only areas with < -200) system(paste("gdal_calc.py -A ", in_path, " --outfile=", out_ras, " --calc=\"(A<=", THRVALUE, ")*A\" --co=\"COMPRESS=LZW\" --type='Int16' --NoDataValue=0 --overwrite", sep=""))
  2. create a binary mask of the raster system(paste("gdal_calc.py -A ", out_ras, " --outfile=", out_mask, " --calc=\"A<0\" --co=\"COMPRESS=LZW\" --type='Byte' --NoDataValue=0 --overwrite", sep=""))
  3. polygonize this mask system(paste("gdal_polygonize.py", out_mask, out_gpkg, sep=" "))
  4. calculate area and filter small polygons
  library(sf)
  diffmask_sf = read_sf(out_gpkg)
  # calculate area
  diffmask_sf$area = diffmask_sf$area <- round(st_area(diffmask_sf))
  #filter out polygons smaller than AREA_THRESHOLD
  diffmask_sf = diffmask_sf[which(diffmask_sf$area > AREA_THRESHOLD),]

Edit: I'm quite happy with a solution with gdal_sieve (see below). However an approach with morpholigcal image processing might even be better, not only allowing to filter by area but also filtering by shape (thin lines), closing holes and connecting polygons that are separated only by thin lines - which suit my use case.

1 Answer 1

1

gdal_sieve (doc) seems to do this efficiently, but doesn't work when polygons are surrounded by NAs (which is the when creating a binary mask with gdal). I circumvented this by manually changing the NA value after creating the binary mask with gdal_calc.

The (added) workflow then is:

  # change no data value, since sieve doesn't filter areas surrounded by NA
  system(paste("gdal_edit.py -a_nodata 255", out_mask, sep=" ")) 
  # sieve with -st AREA_THRESHOLD
  system(paste("gdal_sieve.py -st 5", out_mask, out_mask2, sep=" "))
  # output is always Integer; sieve doesn't allow creation options, set NA value to 0 again
  system(paste("gdal_translate -ot Byte -a_nodata 0 -co \"COMPRESS=LZW\" ", out_mask2, out_mask3, sep=" "))

This is surprisingly fast (a few minutes) and drastically reduces the number of polygons created by the subsequent gdal_polygonize, lowering my computation time from days to minutes.


Update 03-2023: I still believe the solution above with gdal is probably the fastest option. However since I switched to terra for a few other applications/scripts, I am positively surprised how much that sped up raster processing. I wanted to add that, since maybe you are more comfortable looking for a terra solution in R rather than pasting gdal commands to console ;)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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