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)

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. — jdmcbr, Jan 06 '17 at 04:28
Yep - `mask` is a boolean array (as a result of thresholding). — Cate, Jan 06 '17 at 10:08

score 8 · Accepted Answer · answered Jan 06 '17 at 19:51

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)

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). — Cate, Jan 09 '17 at 12:20
I guess my example doesn't integrate the plotting part of your question. You could take a look at this question (http://gis.stackexchange.com/questions/193653/plot-shapefile-on-top-of-raster-using-plot-and-imshow-from-matplotlib) about plotting a shapefile on a raster in python. — jdmcbr, Jan 09 '17 at 17:32

score 2 · Answer 2 · answered Mar 08 '20 at 11:40

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.

answered Mar 08 '20 at 11:40

François

65
6

This approach will create polygons which will not give correct results always because contours are not always polygons and it will always create polygons. – Shubham_geo Dec 21 '20 at 18:48

score 0 · Answer 3 · answered Jan 05 '17 at 17:09

0

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.

answered Jan 05 '17 at 17:09

Jeru Luke

20,118
13
80
87

score 0 · Answer 4 · answered Jan 06 '17 at 19:52

0

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

answered Jan 06 '17 at 19:52

soupault

6,089
4
24
35

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? – Cate Jan 09 '17 at 12:37
@Cate sorry, I'm not familiar with the mentioned formats (shapefile, geojson). I think you should take a look at https://pypi.python.org/pypi/geojson and https://pypi.python.org/pypi/pyshp. Perhaps, the functionality you are looking for is already there. – soupault Jan 09 '17 at 13:01

Shubham_geo · Answer 5 · 2020-12-21T18:50:34.577

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

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

5 Answers5