Polygonize Raster to Shapely Polygons – How to Polygonize Raster to Shapely Polygons

polygonizepython-2.7shapely

I am seeking an open-source python solution to convert raster to polygon (no ArcPy).

I did know the GDAL function to convert raster to polygon, and here is the manual: http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#polygonize-a-raster-band

Nevertheless, I expect that the output can be shapely polygons or any object, temporarily in memory, not saved as a file. Is there any package or code to handle this issue?

If the raster has been processed in a numpy array, the approach is listed below.

Best Answer

Use rasterio of Sean Gillies. It can be easily combined with Fiona (read and write shapefiles) and shapely of the same author.

In the script rasterio_polygonize.py the beginning is

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('a_raster') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))

The result is a generator of GeoJSON features

 geoms = list(results)
 # first feature
 print geoms[0]
 {'geometry': {'type': 'Polygon', 'coordinates': [[(202086.577, 90534.3504440678), (202086.577, 90498.96207), (202121.96537406777, 90498.96207), (202121.96537406777, 90534.3504440678), (202086.577, 90534.3504440678)]]}, 'properties': {'raster_val': 170.52000427246094}}

That you can transform into shapely geometries

from shapely.geometry import shape
print shape(geoms[0]['geometry'])
POLYGON ((202086.577 90534.35044406779, 202086.577 90498.96206999999, 202121.9653740678 90498.96206999999, 202121.9653740678 90534.35044406779, 202086.577 90534.35044406779))

Create geopandas Dataframe and enable easy to use functionalities of spatial join, plotting, save as geojson, ESRI shapefile etc.

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

Related Solutions

Convert Lines to Polygons with Shapely – Guide to Using Polygonize

Create a convex hull from your polyline: https://stackoverflow.com/questions/2964751/how-to-convert-a-geos-multilinestring-to-polygon-using-python

Of course the line has to resemble polygon. If it is eg a river, convex hull will not look very nice.

import fiona
from shapely.geometry import shape, mapping

file = r'e:\...\line.shp'
output = r'e:\temp\output_poly.shp'


schema = {
    'geometry': 'Polygon',
    'properties' : {'id':'int'}
}

with fiona.open(file) as in_file, fiona.open(output, 'w', 'ESRI Shapefile', schema) as out_file:

    for index, row in enumerate(in_file):
        line = shape(row['geometry'])
        hull = line.convex_hull
        out_file.write({
            'geometry': mapping(hull),
            'properties': {'id': index},
        })
        print(hull)

Also if you want to do it manually, you can. You just have to create polygon from vertexes and add one vertex, which is the same as the first one from line so you have a closed ring.

edit: Since convex hull is not working for you, try this piece of code. It will take the points from Line and create polygon, it will put the first point as last point of line so it is closed ring.

import fiona
from shapely.geometry import shape, mapping, LineString, Polygon

file = r'e:\Temp\lines.shp'
output = r'e:\temp\output.shp'


schema = {
    'geometry': 'Polygon',
    'properties' : {'id':'int'}
}

with fiona.open(file) as in_file, fiona.open(output, 'w', 'ESRI Shapefile', schema) as out_file:
    for index_line, row in enumerate(in_file):
        line = shape(row['geometry'])
        coordinates = []

        if isinstance(line, LineString):
            for index, point in enumerate(line.coords):
                if index == 0:
                    first_pt = point
                coordinates.append(point)

            coordinates.append(first_pt)
            if len(coordinates) >= 3:
                polygon = Polygon(coordinates)
                print(polygon)
                out_file.write({
                    'geometry': mapping(polygon),
                    'properties': {'id': index_line},
                })

ArcGIS – Converting Dataframe with WKT Column to Spatial Dataframe in Python API

Looks like this issue was reported to the arcgis-python repo and has been fixed but not yet released.

As a workaround you may consider rolling the ArcGIS Python package back to 1.9.1 rather than the current 2.0.0 release.

You may also try setting arcgis.geometry._types._HASSHAPELY = True.

Best Answer

Related Solutions

Convert Lines to Polygons with Shapely – Guide to Using Polygonize

ArcGIS – Converting Dataframe with WKT Column to Spatial Dataframe in Python API

Related Question