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?

  • What's wrong with raster? Are you sure that drawing a bunch of polygons is going to be faster than rendering a raster? If you're keen then you need to use the OGR component of GDAL to access the results from rasterize, have a look at the code, it's online and in C++ searchcode.com/codesearch/view/13193937 – Michael Stimson Dec 11 '18 at 7:53
  • @MichaelStimson, I'm using a rendering engine that is is optimized for polygons and not for rendering raster – Graviton Dec 11 '18 at 8:34
  • Have you looked at the OGR tutorial? gdal.org/ogr_apitut.html GDAL_Polygonize creates a layer (probably a shapefile) which can be opened by OGR then iterated with layer->GetNextFeature() (for values and geometry) then iterate the geometry obtained with OGRfeature->GetGeometryRef->toPoint as OGRPoint* which have X and Y to make up your list. – Michael Stimson Dec 12 '18 at 6:36
  • @MichaelStimson, So let me understand this: one OGRFeature refers to one polygon? Because this is the impression I get when I play with the API... – Graviton Dec 12 '18 at 8:59
  • 1
    One OGRFeature is one row in the feature class which contains one geometry as well as all the associated attributes. Because a feature class (layer) usually contains more than one row (feature) it is necessary to cursor through them; beware though rows are returned in no particular order, usually OID ascending but don't build code that relies on that. When GetNextFeature() returns NULL you have reached the end of the table (layer). – Michael Stimson Dec 13 '18 at 2:18

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

  • I have remarked in the question that I have seen the cookbook reference , and why it is of no help to me( or rather , I don't know how to make use of it), and that's why I ask the question here – Graviton Dec 12 '18 at 1:44
  • Sorry, I guess I'm missing something here. Are you rendering each individual feature using a custom made code? In other words, is the question about polygonizing or how to access polygons? – Alessio Arena Dec 12 '18 at 4:12
  • Yes, I'm rendering using my own custom made code. – Graviton Dec 12 '18 at 4:20
  • This is a question about how to render the raster more efficiently; I assume that by grouping many raster cells into a few large polygons will make the whole rendering process faster, in my own rendering engine. So I am asking how to actually group the adjacent raster cells with the same values into one polygon – Graviton Dec 12 '18 at 4:21
  • I modified the answer pointing to additional examples on how to access each individual polygons. Regarding the grouping, that is actually already done by Polygonize. Practically you should be able to run polygonize, then iterate through the layer and get your already grouped features. – Alessio Arena Dec 12 '18 at 4:30

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.