Python分类遥感图像并提取矢量边界

图像来源于Python裁剪遥感图像中得到的裁剪图像

# 参考《python地理空间分析指南》
import numpy as np
import cv2
from PIL import Image
from osgeo import gdal_array, gdal, ogr, osr
import shapefile
import pngcanvas
import matplotlib.pyplot as plt

# 先使用分类的方法将图像分为两类
src = 'phase2/maskcliped.tif'  # 分类后的原图
tgt = 'phase2/class_2.tif'  # 保存图像

src_arr = gdal_array.LoadFile(src)
# 根据类别数将直方图分割成2个颜色区间
classes = gdal_array.numpy.histogram(src_arr, bins=2)[1]
# 颜色查表的记录数必须为len(classes)+1
lut = [[255,0,0], [255,255,255], [0,0,0]]
# 分类初始值
start = 1
# 创建RGB颜色的JPEG输出图片
rgb = gdal_array.numpy.zeros((3, src_arr.shape[1], src_arr.shape[2]), gdal_array.numpy.float32)
# 处理所有类并申明颜色
for i in range(len(classes)):
    mask = gdal_array.numpy.logical_and(start <= src_arr, src_arr <= classes[i])  # 这里是3通道的
    masked = ((mask[0] == mask[1]) == mask[2])  # 3通道bool型合并为1通道
    for j in range(len(lut[i])):
        rgb[j] = gdal_array.numpy.choose(masked, (rgb[j], lut[i][j]))  # 需要1通道的
    start = classes[i] + 1
print(rgb.shape)
# 保存图像
out = gdal_array.SaveArray(rgb.astype(gdal_array.numpy.uint8), tgt, format='GTIFF', prototype=src)
out = None

(3, 986, 951)

# 加载原图及分类结果
NUMS = 65536
raw_arr = cv2.merge((src_arr[0]/float(NUMS), src_arr[1]/float(NUMS), src_arr[2]/float(NUMS)))
cls_img = Image.open(tgt)
plt.figure(figsize=(20, 10))
plt.subplot(121),plt.imshow(raw_arr)
plt.subplot(122),plt.imshow(cls_img)
plt.show()

# 阈值化输入的栅格图像
src_c = 'phase2/class_2.tif'
tgt_c = 'phase2/extract.shp'  # 输出的shp
tgt_layer = 'extract'  # 图层名称
# 打开栅格图像
src_ds = gdal.Open(src_c, gdal.GA_ReadOnly)
# 读取第一波段
band_1 = src_ds.GetRasterBand(1)
# 使用该层作为遮罩
shp_mask = band_1.GetMaskBand()
# 创建shp文件
driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.CreateDataSource(tgt_c)
# 拷贝空间索引
srs = osr.SpatialReference()
srs.ImportFromWkt(src_ds.GetProjectionRef())
layer = shp.CreateLayer(tgt_layer, srs=srs)
# 创建dbf
fd = ogr.FieldDefn('DN', ogr.OFTInteger)
layer.CreateField(fd)
dst_field = 0
# 从图片提取特征
extract = gdal.Polygonize(band_1, shp_mask, layer, dst_field, [], None)
extract = None
shp = None  # 一定要记得关闭（之前没关闭shp一直打不开是空白）
print('ok')

ok

# 显示输出结果
png_path = 'phase2/extract.png'
r = shapefile.Reader(tgt_c)
x_dist = r.bbox[2] - r.bbox[0]
y_dist = r.bbox[3] - r.bbox[1]
iwid = 951
ihei = 986
x_rat = iwid / x_dist
y_rat = ihei / y_dist
polygons = []
for shape in r.shapes():
    for i in range(len(shape.parts)):
        pixels = []
        pt = None
        if i < len(shape.parts) - 1:
            pt = shape.points[shape.parts[i]:shape.parts[i+1]]
        else:
            pt = shape.points[shape.parts[i]:]
        for x,y in pt:
            px = int(iwid - (r.bbox[2] - x) * x_rat)
            py = int((r.bbox[3] - y) * y_rat)
            pixels.append([px, py])
        polygons.append(pixels)
c = pngcanvas.PNGCanvas(iwid, ihei)
for p in polygons:
    c.polyline(p)
f = open(png_path, 'wb')
f.write(c.dump())
f.close()

img = Image.open(png_path)
plt.figure(figsize=(10, 10))
plt.imshow(img)
plt.show()

用ArcGIS打开的效果如下