How to polygonize raster to shapely polygons
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)
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
#segimg=glob.glob('Poly.tif')[0]
#src_ds = gdal.Open(segimg, GA_ReadOnly )
#srcband=src_ds.GetRasterBand(1)
#myarray=srcband.ReadAsArray()
#these lines use gdal to import an image. 'myarray' can be any numpy array
mypoly=[]
for vec in rasterio.features.shapes(myarray):
mypoly.append(shape(vec))
The way to install rasterio is by 'conda install -c https://conda.anaconda.org/ioos rasterio', if there is an installation issue.