Constructing polygons from GDAL Polygonize?

Question

I have a geotiff file, which, for each cell ( pixel), it contains a discrete value ( you can think of this value as the "fitness" value of the cell). For the neighboring cells, they can consist of the same discrete value, but not always.

I want to create a polygon for all the adjacent cells that share the same discrete value. A polygon is defined by the following usual class ( in pseudo code)

class Point
{
    double X;
    double Y;
}

class Polygon
{
    List<Point> Pts;
}

The reason for this is that, for my rendering engine, I would rather render a few large polygons, rather than all the individual raster cells one by one for performance consideration.

I believe that GDALPolygonize is the way to go. However I don't quite know how to turn the output from the function into a list of Polygons that I can use.

The example from the cookbook straightaway makes use of a library to render the output file, which is not suitable for my purpose.

I mainly working with GDAL ( C++) or GDAL.Net (C#), but examples in other languages are acceptable as the syntax is quite close with C++ version.

Or is there other, more direct method that I can use, both in or out of GDAL?

Answer 1

Had a look at this link? It is the python gdal cookbook, but can give you some close references to your use case.

This operation will write your polygons into a layer. In python you can access those individual polygons and features using those examples. If you look at the GDAL documentation you will see that the object returned by GDALPolygonize belongs to the class OGRLayer, and that has similar methods in C (GetFeatureCount and GetFeature)

Just to bring some clarity around what you are accessing and when: In this cookbook example you are initialising a new shapefile, that can be commonly referred to as a FeatureCollection. Here is how they do it

drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )

A feature collection may have one or more FeatureLayers. In this example you initialise this by doing:

dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

At this point the layer is empty, but that is the destination of your Polygonize output. This operation will populate the layer with new Features (in your case polygons).

For your application you want to access this FeatureLayer and iterate through each Feature to access its details and characteristics