I am working on a workflow to identify tree stand boundaries. My desired result is the blue lines - I know it will never be like that, but I need it to be as good as possible. My best result is the red lines. The final result must not be polygons, but lines.

enter image description here

So this is what I have done so far: I am using a Normalized DSM with pixel-size 2x2 as my starting layer. This is my workflow: enter image description here

I use i.segment where my best results have been with Difference Threshold=0.4 and Minimum Number of cells = 150. Then I use r.to.vect because Polygonize for some reason didn't work in this workflow (?). After this I convert to lines, and delete duplicates. I use line simplification with Tolerence 8.5 (best result) And last I use v.generalize with Maximal Tolerance Value=21 with the Douglas algorithm.

So I use four variables to make my result better. I have tried to use some more of the other variables, but honestly I am not sure how they affect the result.

So I am asking here for any suggestions to make my result better?

Update: I managed to make it slightly better by keep changing the variables above. I tried some other variables:

  • In segmentation I changed the amount of memory to use in MB, but it had no effect.
  • In v.generalize i tried all different algorithms, and Douglas appears to be the best.
  • I also tried changing the Look-ahead parameter in v.generalize, but it had no positive effect.

So I am still hoping for any "out-of-the-box" suggestions to make it better.

Update-2: I am still working on this. I have tried different variations for raster inputs. I have tried the following:

  1. Normalized DSM (resolution 2x2)
  2. Normalized DSM (2x2) + orthofoto
  3. Normalized DSM (2x2) + orthofoto + Orthofoto CIR
  4. Orthofoto + Orthofoto CIR
  5. Normalized DSM (2x2) + DSM (2x2)

I have also tried different resolutions, but 2x2 gives the best result. I have also learned that combining three raster inputs gives a worse result. My best results is either just using Normalized DSM (2x2) or the combined input with Normalized DSM (2x2) + DSM (2x2).

Still hoping for some suggestions. I am open for trying an other open source product.

  • Right now I am trying with more inputs,as you wrote. But its back to zero in many ways. So I am working with that. If you have an other open source solution, I am open to try it. Sep 16, 2020 at 13:54

1 Answer 1


I had a go with some simple Python. Hope this is still acceptable. Basically, I did a bit of blurring (Gaussian smoothing, you can blur the image more or less to make the homogeneous areas more homogeneous), and then applied a segmentation algorithm. My script is heavily based on some scikit-image documentation, and this is what I get in terms of segments:

segmentation results

You can of course tweak the parameters, and if you wanted to save the output as an e.g. GeoTIFF you could then extract polygons from, just save the segments_slic as GeoTIFF with the appropriate geolocation. I left that bit of the code commented out at the end.

# Big import block
from PIL import Image
import matplotlib.pyplot as plt
import numpy as np

from skimage.color import rgb2gray
from skimage.segmentation import slic
from skimage.segmentation import mark_boundaries

# Open image
im_fp = Image.open("forest.png")
# Turn into numpy array, discard last band
# (alpha channel?)
img = np.asarray(im_fp)[:, :, :-1]

# Smooth spatially a bit
img=ndi.gaussian_filter(img, sigma=5)

# Calculate segments. 
# Feel free to tweak n_segments, compactness
# and sigma
segments_slic = slic(img, n_segments=25, compactness=50, 
                     sigma=1, start_label=1)

# Plot image and boundaries
plt.imshow(mark_boundaries(img, segments_slic))

# Export as GeoTIFF. Check
#drv = gdal.GetDriverByName('GTiff')     # create driver for writing geotiff file
#outRaster = drv.CreateCopy('forest_segments.tif', 'forest.tif, , 0 )   # create new copy of inut raster on disk
#outRaster.SetGeoTransform([top_left_x, hor_resolution, 0, top_left_y, 0, -vert_res])
#newBand = outRaster.GetRasterBand(1)                               # get the first (and only) band of the new copy
#newBand.WriteArray(segments_slic)                                           # write array data to this band 
#outRaster = None

  • 1
    any rule of thumb on how to optimise the choice of SLIC parameters (e.g. segments, compactness)?
    – Nico
    Feb 15 at 18:29
  • 2
    n_segments is roughly the number of patches on the image, compactness is basically a tradeoff between square pixels and the underlying data. Choose a lower compactness for less blocky patches, and viceversa. sigma: smoothing kernel width. I already applied a sigma of 5, but maybe not needed? All options here
    – Jose
    Feb 15 at 19:01

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.