I am facing with the problem of having a set of .grc files having just one band. However, this file contains several categories that belong, in the end, to different portion of the image (I guess).
When I generate the shape file out of this (or directly a PostGIS geometry), all these categories get merged into a single file having the separated regions but DN 1, DN 2, DN 3 .. as names for them. I would like to preserve the category names in the shape file.
This command generates currently the shape files:
gdal_polygonize.py coverage.grc -f "ESRI Shapefile" output.shp
Are my assumptions correct? How can I do that?
I managed to solve this by creating a custom script in Python:
catnames = rasterband.GetCategoryNames()
memdriver = ogr.GetDriverByName( 'Memory' )
dst_ds = memdriver.CreateDataSource( "out" )
# create a layer with a DN field that will contain the index of
# each category automatically
dst_layer = dst_ds.CreateLayer("out", srs = None )
fd = ogr.FieldDefn( 'DN', ogr.OFTInteger )
dst_layer.CreateField(fd)
gdal.Polygonize( rasterband, None, dst_layer, 0)
# extract features map
featmap = {}
for feature in dst_layer:
# use the field defined above, create a mapping between
# index -> feature name
featname = feature.GetField("DN")
if featname not in featmap.keys():
featmap[featname] = [feature]
else:
featmap[featname].append(feature)
From this point on, dst_layer contains the vector shape that can be easily saved to shape using:
driver = ogr.GetDriverByName("ESRI Shapefile")
shapeData = driver.CreateDataSource(args.output_folder) #so there we will store our data
for featidx, features in featmap.iteritems():
featname = catnames[featidx]
layer = shapeData.CreateLayer(featname, spatialReference, ogr.wkbPolygon)
for feature in features:
layer.CreateFeature(feature)
shapeData.Destroy() #lets close the shapefile
This will create a .shp file for each feature name (category). Of course, this can be tuned to create one single file but this is out of the scope of my question.