如何将在 scikit-image find_contours 中创建的轮廓导出到 shapefile 或 geojson？

Question

我正在尝试在卫星图像运行之后将 scikit-image.measure.find_contours() 函数的结果导出为 shapefile 或 geojson。

输出是一个类似（行，列）的数组，坐标沿等高线，其中有很多。

如何绘制各种等高线的坐标，并将其导出到 shapefile（可以设置适当的投影等）？

我当前的代码，其中 'mask' 是我处理过的图像：

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

您应该安装 python 库 geojson 并使用它。

为了使用图像中的坐标和标记对象，您应该使用库 shapely。

Answer 2

根据 rasterio 和 fiona 的主要开发人员 post 改编的以下内容应该可以工作，但我相信您需要再适应一点。它使用 rasterio.features.shapes 来识别图像中具有某些值的连续区域和 return 相关联的坐标，基于栅格的变换。然后使用 fiona.

将这些记录写入 shapefile

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 3

@Cate 你可以使用那些 row, column 坐标矩阵并通过 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

简单地使用skimage:

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

i 是您要绘制覆盖尺寸 image_shape 原始图像的轮廓的索引。

Answer 5

这是我的食谱，效果很好。

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

如何将在 scikit-image find_contours 中创建的轮廓导出到 shapefile 或 geojson？

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

python

arrays

shapefile

contour

scikit-image