How to export contours created in scikit-image find_contours to shapefile or geojson?

Question

I'm trying to export the results of the scikit-image.measure.find_contours() function as a shapefile or geojson after running on a satellite image.

The output is an array like (row, column) with coordinates along the contours, of which there are many.

How do I plot the coordinates of the various contours, and export this to a shapefile (can set appropriate projection etc.)?

My current code where 'mask' is my processed image:

from skimage import measure
import matplotlib.pyplot as plt

contours = measure.find_contours(mask, 0.5)

plt.imshow(mask)
for n, contour in enumerate(contours):
    plt.plot(contour[:,1], contour[:, 0], linewidth=1)

Answer 1

Something along the lines of the following, adapted from a post by the primary developer of rasterio and fiona, should work, though I'm sure you'll need to adapt a little more. It uses rasterio.features.shapes to identify contiguous regions in an image that have some value and return the associated coordinates, based on the transform of the raster. It then writes those records to a shapefile using fiona.

import fiona
import rasterio.features

schema = {"geometry": "Polygon", "properties": {"value": "int"}}

with rasterio.open(raster_filename) as raster:
    image = raster.read()
    # use your function to generate mask
    mask = your_thresholding_function(image)
    # and convert to uint8 for rasterio.features.shapes
    mask = mask.astype('uint8')
    shapes = rasterio.features.shapes(mask, transform=raster.transform)
    # select the records from shapes where the value is 1,
    # or where the mask was True
    records = [{"geometry": geometry, "properties": {"value": value}}
               for (geometry, value) in shapes if value == 1]
    with fiona.open(shape_filename, "w", "ESRI Shapefile",
                    crs=raster.crs.data, schema=schema) as out_file:
        out_file.writerecords(records)

---

Answer 2

You should install python library geojson and work with that.

In order to play with the coordinates and mark objects present in the image you should use the library shapely.

Answer 3

@Cate you could use those row, column coordinate matrices and plot them via http://scikit-image.org/docs/dev/api/skimage.draw.html#skimage.draw.polygon (filled polygon), http://scikit-image.org/docs/dev/api/skimage.draw.html#skimage.draw.polygon_perimeter (only perimeter), or create your custom polygon plotting function on top of http://matplotlib.org/api/patches_api.html#matplotlib.patches.Polygon.

Answer 4

Adding one more solution, for converting the output of skimage.measure.find_contours to a shapefile using pyshp:

import shapefile
import skimage.measure

contours = skimage.measure.find_contours(mask, 0.5)
points = []
parts = []
for loop in contours:
    parts.append(len(points))
    points.extend(geo_loop)
shape = shapefile.Shape(
    shapeType=shapefile.POLYGON, points=points, parts=parts)

with shapefile.Writer(output_shapefile,
                      shapeType=shapefile.POLYGON) as w:
    w.field('Field', 'C')
    w.record('polygon_id')
    w.shape(shape)

Answer 5

This is my recipe and it works quite well.

import skimage
import gdal
import matplotlib.pyplot as plt
import numpy as np
import rasterio
import shapely
import fiona

#Open raster with gdal
image=gdal.Open('A.tif')
im=image.ReadAsArray()

#out variable stores the contours
out=skimage.measure.find_contours(im,0.5)
# Here,0.5 is taken  assuming it is a binary raster
# but the default value is taken as (np.max(im)+np.min(im))/2

fig, ax = plt.subplots()
ax.imshow(im, cmap=plt.cm.gray)
#cs list will contain all the 2D Line objects
cs=[]
for contour in out:
    cs.append(ax.plot(contour[:, 1], contour[:, 0], linewidth=2))

ax.axis('image')
#Show image with contours
plt.show()

#Read band 1 of raster or as per the usage it can be tweaked
with rasterio.open('A.tif') as raster:
    image = raster.read()[0,:,:]
    
#Create list poly containing all the linestrings of contours
from shapely.geometry import mapping,MultiLineString,LineString
poly=[]
for i in cs:
    
    x=i[0].get_xdata()
    y=i[0].get_ydata()
    aa=rasterio.transform.xy(raster.transform,y,x)
    poly.append(LineString([(i[0], i[1]) for i in zip(aa[0],aa[1])]))

#Create a list of wkt strings
list_lstrings =  [shapely.wkt.loads(p.wkt) for p in poly]
# Create a MultiLineString object from the list
mult=shapely.geometry.MultiLineString(list_lstrings)

#Inputting projection info
from fiona.crs import from_epsg
crs = from_epsg(4326)

#Create schema    
schema = {
    'geometry': 'MultiLineString',
    'properties': {'id': 'int'},
}

# Write a new Shapefile
with fiona.open('U:\\new_shape.shp', 'w', 'ESRI Shapefile', schema,crs=crs) as c:
    
    c.write({
        'geometry': mapping(mult),
        'properties': {'id': 1},
    })
    

Answer 6

simply use skimage:

from skimage.draw import polygon2mask mask = polygon2mask(image_shape, contours[i])

i is the index of the contour you want to plot overlaying your original image of dimension image_shape.