How to convert float raster to vector with python GDAL

Question

I have float raster and now I want to convert it to vector. How is it possible with the Python GDAL library?

I have tried with gdal_polygonize.py of GDAL utilities on the command line and it worked excellently. But this utility is based on GDALPolygonize() of C++ library and I want to C++ method GDALFPolygonize() to be used instead, which manipulates float raster data as far as I know.

Answer 1

Try using rasterio, which uses GDALFPolygonize on float arrays.

import numpy as np
import rasterio
from affine import Affine
from shapely.geometry import shape

# triangular array    
ar = np.tri(5, dtype='f')
print(ar)

for shp, val in rasterio.features.shapes(ar, transform=Affine(1, 0, 0, 0, -1, 5)):
    print('%s: %s' % (val, shape(shp)))

shows:

[[ 1.  0.  0.  0.  0.]
 [ 1.  1.  0.  0.  0.]
 [ 1.  1.  1.  0.  0.]
 [ 1.  1.  1.  1.  0.]
 [ 1.  1.  1.  1.  1.]]
1.0: POLYGON ((0 0, 0 5, 5 5, 5 4, 4 4, 4 3, 3 3, 3 2, 2 2, 2 1, 1 1, 1 0, 0 0))
0.0: POLYGON ((1 0, 1 1, 2 1, 2 2, 3 2, 3 3, 4 3, 4 4, 5 4, 5 0, 1 0))

Or visualised with blue for 1.0 and red for 0.0:

Answer 2

You can't use GDALFPolygonize with the GDAL python bindings without modifying the source code and recompiling as it isn't exposed in the GDAL swig interface.

Note: as at Feb 2016, GDALFPolygonize IS exposed in the GDAL SVN trunk source, but is not in either of the latest releases (1.11.4/2.0.2).

To polygonize your raster, you will need to convert from float to integer. If you want to retain some decimal places multiply your raster by 10^N where N is the number of decimal places you want to retain. For example, to keep 3 decimal places multiply by 10^3 = 1000.

# Multiply raster values by 1000
gdal_calc.py -A raster.tif --outfile=raster1000.tif --calc="A*1000"

# Convert to vector
gdal_polygonize.py raster1000.tif vector1000.shp

If you want to convert your polygon attributes back to float then just divide by the same value.

# Divide vector values by 1000
ogr2ogr vector.shp vector1000.shp -sql "select (cast(DN as float) * 0.001) as DN FROM vector1000"

Note there is no point converting floating point rasters which represent continuous surfaces to polygons as you will get pretty much 1 polygon per pixel which is very inefficient. Almost any analysis you can think of with such data is much more efficient if you leave the data in raster format.

Answer 3

A bit late, but commenting if anyone finds this useful. Let's say you have a float raster (raster_path) and you want to convert it to shapefile. You can use gdal FPolygonize within python as:

source_raster = gdal.Open(raster_path)      
band = gdal.Open(raster_path).GetRasterBand(1)
driver = ogr.GetDriverByName('ESRI Shapefile')        
out_data = driver.CreateDataSource(shapefile_path)
# getting projection from source raster
srs = osr.SpatialReference()
srs.ImportFromWkt(source_raster.GetProjectionRef())
# create layer with projection
out_layer = out_data.CreateLayer(raster_path.split('.')[0], srs)        
new_field = ogr.FieldDefn('field_name', ogr.OFTReal)
out_layer.CreateField(new_field)        
gdal.FPolygonize(band, band, out_layer, 0, [], callback=None)        
out_data.Destroy()
source_raster = None