[SOLVED] How to export contours created in scikit-image find

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

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)

Solution

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)