5

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.

5

Try using rasterio, which uses GDALFPolygonize on float arrays.

import numpy as np
import rasterio.features
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:

from JTS

5

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.

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.