[GIS] How to convert float raster to vector with python GDAL

gdalpythonvectorization

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.

Best Answer

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:

from JTS

Related Solutions

[GIS] Creating vector from raster using GDAL

As Matthew asked and then I commented on, I apparently was confused about the world file creation and in which step to do that in.

Adding -co WORLDFILE=YES in step 3 above created a .wld file for the resulting PNG image. I then used gdal_contour, brought into QGIS, and it worked!

[GIS] Convert points to lines with Python GDAL

Why not work globally ?

calculate the distances between all points
union the resulting lines pointx - pointy with a distance < 14m

I will use Shapely, much easier for resolving these kinds of problems. You must iterate through all pairs of points to calculate the distance once (as distance point1-point2 = distance point2-point1). There are many solutions in Python and I choose the itertools standard module with combinations.

example:

myPointDict = {0:(1,1), 1:(2,2), 2:(3,3),3:(4,4),4:(5,5)}
import itertools
for i in  itertools.combinations(PointDict.values(), 2):
   print i
((1, 1), (2, 2))
((1, 1), (3, 3))
((1, 1), (4, 4))
((1, 1), (5, 5))
((2, 2), (3, 3))
((2, 2), (4, 4))
((2, 2), (5, 5))
((3, 3), (4, 4))
((3, 3), (5, 5))
((4, 4), (5, 5))

With ogr (look at the Python GDAL/OGR Cookbook!):

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(x,y)
distance =  point1.Distance(point2)
line = ogr.Geometry(ogr.wkbLineString) 
    line.AddPoint(x1, y1)
    ....
    line.AddPoint(xn,yn)

With shapely:

point = Point(x,y)
distance = Point(x1,y1).distance(Point(x2,y2)
linestring = LineString([point1,..., pointn]

So, in your case:

from shapely.geometry import Point, LineString
# creation of a empty line for unioning the resulting geometries
line = LineString()
for i in  itertools.combinations(pointDict.values(), 2):
     # if distance < 14m union the line ptx-pty to line
     if Point(i[0]).distance(Point(i[1])) < 14:
            line = line.union(LineString([(Point(i[0]).x, Point(i[0]).y), (Point(i[1]).x, Point(i[1]).y)]))
     # result 
     print line.wkt
     'MULTILINESTRING ((345672.493225679441821 1267286.555012494325638,345681.57590266619809 1267286.555012494325638),(345672.493225679441821 1267286.555012494325638,345663.410548692685552 1267277.472335507394746),(345672.493225679441821 1267286.555012494325638,345681.57590266619809 1267277.472335507394746),(345681.57590266619809 1267286.555012494325638,345681.57590266619809 1267277.472335507394746),(345654.327871705929283 1267277.472335507394746,345663.410548692685552 1267277.472335507394746),(345654.327871705929283 1267277.472335507394746,345645.245194719173014 1267268.389658520696685),(345681.57590266619809 1267277.472335507394746,345690.658579652954359 1267268.389658520696685),(345636.162517732358538 1267268.389658520696685,345645.245194719173014 1267268.389658520696685),(345636.162517732358538 1267268.389658520696685,345627.079840745602269 1267259.306981533998623),(345690.658579652954359 1267268.389658520696685,345681.57590266619809 1267259.306981533998623),(345608.914486772089731 1267259.306981533998623,345617.997163758846 1267259.306981533998623),(345608.914486772089731 1267259.306981533998623,345599.831809785333462 1267250.224304547300562),(345617.997163758846 1267259.306981533998623,345627.079840745602269 1267259.306981533998623),(345681.57590266619809 1267259.306981533998623,345672.493225679441821 1267250.224304547300562),(345590.749132798577193 1267250.224304547300562,345599.831809785333462 1267250.224304547300562),(345672.493225679441821 1267250.224304547300562,345663.410548692685552 1267241.14162756036967))'

And if you want to use the end of your script:

multiline = ogr.CreateGeometryFromWkt(line.wkt)

or using Fiona (an easier Python wrapper of the ogr library)

import fiona
from shapely.geometry import mapping
# schema of the shapefile
schema = {'geometry': 'MultiLineString','properties': {'test': 'int'}}
with fiona.open('myshp3.shp','w','ESRI Shapefile', schema) as c:
       record = {'geometry':mapping(line), 'properties':{'test':1}}
       c.write(record)

Result:

enter image description here

But, I do not know if this is what you want.

Best Answer

Related Solutions

[GIS] Creating vector from raster using GDAL

[GIS] Convert points to lines with Python GDAL

Related Question