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.

  • I am not sure if one mix those two packages. Are you aware of mutual incompatibilites between rasterio and gdal? > Mutual Incompatibilities Rasterio and GDAL’s bindings can contend for global GDAL objects. Unless you have deep knowledge about both packages, choose exactly one of import osgeo.gdal or import rasterio see here – Nerolf05 Mar 15 at 13:53

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)
  • if the raster has been processed as a numpy array, is there any way to convert a numpy array as polygons? Thanks! – Vicky Liau Apr 3 '16 at 18:45
  • In theory, yes- – gene Apr 3 '16 at 19:17
  • 1
    The mask variable and parameter in your example seems unnecessary. I would however recommend to add if value > src.nodata to the list comprehension to make use of the source's nodata value and discard any shapes that correspond to it. Not sure what would happen if there is no nodata vaue though. :o) – bugmenot123 Apr 28 '17 at 18:58
  • 6
    meanwhile they changed rasterio.drivers to rasterio.Env and src.affine to src.transform – Leo Dec 3 '18 at 9:31
  • When I tried this I get a polygon with a crs of None. So I tried setting the crs to the original image's but then it shows it in the completely wrong place. Any ideas? – Emtomp Apr 1 '20 at 14:50

Here is my implementation.

from osgeo import ogr, gdal, osr
from osgeo.gdalnumeric import *  
from osgeo.gdalconst import * 
import fiona
from shapely.geometry import shape
import rasterio.features

#src_ds = gdal.Open(segimg, GA_ReadOnly )
#these lines use gdal to import an image. 'myarray' can be any numpy array

for vec in rasterio.features.shapes(myarray):

The way to install rasterio is by 'conda install -c https://conda.anaconda.org/ioos rasterio', if there is an installation issue.

  • The result of rasterio is directly a numpy array, therefore you don't need myarray=srcband.ReadAsArray() #or any array – gene Apr 4 '16 at 6:16
  • @gene I revised the note. This line (myarray=srcband.ReadAsArray()) uses gdal to import the image. – Vicky Liau Apr 4 '16 at 18:21
  • import the image as a numpy array and rasterio import directly the image as a numpy array – gene Apr 4 '16 at 18:51
  • This worked for me although I had to index vec because it was being returned as a tuple. shape(vec[0]) – narmaps Oct 9 '19 at 17:00

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.