我试图将scikit-image.measure.find_contours()函数的输出作为一个shapefile或geojson导出，该函数在卫星图像上运行后产生(row,column)格式的坐标数组，表示轮廓线上的坐标点。
如何绘制各个轮廓线的坐标，并将其导出为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)
以下翻译参考自rasterio和fiona的主要开发者文章，应该可以满足需求，不过你可能需要做一些修改。它使用rasterio.features.shapes在图像中识别具有某些值的连续区域，并基于栅格变换返回相关坐标。然后使用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)
rasterio.open()代替我原来的gdal.Open()是可以的，因为我可以直接对数组进行任何图像操作（阈值化和其他几个步骤）。在rasterio.features.shapes之前添加mask = np.invert(mask)允许我提取反向区域（在这种情况下对我很有用）。 - Cate只需使用skimage：
from skimage.draw import polygon2mask
mask = polygon2mask(image_shape, contours[i])
i是你想要在具有image_shape维度的原始图像上叠加绘制的轮廓的索引。
你需要安装Python库geojson并使用它。
为了处理坐标和标记图像中的对象，你应该使用shapely库。
@Cate 你可以使用那些 行、列 坐标矩阵，并通过 http://scikit-image.org/docs/dev/api/skimage.draw.html#skimage.draw.polygon（填充多边形）、http://scikit-image.org/docs/dev/api/skimage.draw.html#skimage.draw.polygon_perimeter（仅绘制轮廓）或在 http://matplotlib.org/api/patches_api.html#matplotlib.patches.Polygon 上创建自定义多边形绘图函数。
find_contours() 输出的嵌套列表结构中使其工作。对于每个多边形子列表，您会创建一个 shapefile 并在迭代列表时将其附加到其中，还是有一种简单的方法可以一次性将它们全部写入？ - Cate这是我的配方，而且它的效果相当不错。
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},
})
mask是一个布尔数组吗？如果是的话，可能有一种更直接的方法可以得到您要寻找的结果。 - jdmcbrmask是一个布尔数组（通过阈值处理得到）。 - Cate