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)

Deforest

Given the question you asked yesterday, is mask a boolean array? If so, there might be a slightly more straightforward way of getting to the result you're looking for. – Sesame

Yep - mask is a boolean array (as a result of thresholding). – Deforest

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)

---

Sesame

This works with very few changes. Using rasterio.open() instead of my original gdal.Open() is fine as I can go straight into any image manipulation on the array (thresholding and a few other steps). Adding mask = np.invert(mask) just before rasterio.features.shapes allows me to pull out the inverse regions (useful for me in this case). – Deforest

I guess my example doesn't integrate the plotting part of your question. You could take a look at this question (gis.stackexchange.com/questions/193653/…) about plotting a shapefile on a raster in python. – Sesame

Answer 2

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.

Gally

This approach will create polygons which will not give correct results always because contours are not always polygons and it will always create polygons. – Output

Answer 3

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.

Laciniate

Answer 4

@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.

Karat

Thanks @soupault, trying to get this working within the nested list structure that find_contours() puts out. With each polygon as a sublist would you create a shapefile and append to it as you iterate through the list, or is there an easy way to write them all at once? – Deforest

@Deforest sorry, I'm not familiar with the mentioned formats (shapefile, geojson). I think you should take a look at pypi.python.org/pypi/geojson and pypi.python.org/pypi/pyshp. Perhaps, the functionality you are looking for is already there. – Karat

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},
    })
    

Output

Answer 6

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)

Aston