skimage进行超像素分割创建兴趣区的方法

前段时间写了一篇遥感影像面向对象分类的python实现，意外收到了读者发来的邮件，问我能不能直接按照选定的研究区进行分割，办法肯定比困难多，其实我也早就有想法做这个的实现，方法很简单，在skimage包里就有集成，所以这里简单介绍一下。

方法实现

主要是基于from skimage import morphology方法，morphology里的remove_small_holes函数。

函数用法：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37



from skimage.segmentation import slic
from skimage import morphology

def create_mask(input_arra):
    mask = morphology.remove_small_holes(morphology.remove_small_objects(input_arra > 0, 500),500)  
    #这里算是mask的关键，tif影像中通常以0划分有效区和无效区，因此这里设置‘input_arra’参数为>0，其他基本就是通用写法。

    mask = morphology.opening(mask, morphology.disk(3))
    return mask


def segementation_img(input_raster, output_raster):
  
    dataset = gdal.Open(input_raster)
    im_width = dataset.RasterXSize
    im_height = dataset.RasterYSize
    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height) #0, 0, im_width, im_height
    bandnum = dataset.RasterCount
    im_data = im_data.astype(int)
    if bandnum == 1:
        temp, mask_arra = im_data.transpose((1,0))
    
    elif bandnum == 2:
        im_data = im_data[0]
        temp, mask_arra = im_data.transpose((1,0))
    else:
        im_data = im_data[0:3]
        temp = im_data.transpose((2, 1, 0))
        mask_arra = temp[:, : , 0]
    mask = create_mask(mask_arra)
    seg_func = slic(temp, n_segments=2000, compactness=10, mask=mask)  #1
    #seg_func = quickshift(temp, ratio=1.0, kernel_size=5)
    label = seg_func.transpose((1,0))
    write_img(output_raster, im_proj, im_geotrans, label)

有了以上函数，那怎么利用输入的shp进行区域分割呢？这里还涉及了shp转栅格并二值化，有兴趣的朋友还可以继续挖掘这一步的写法😙，提示一下，基于gdal，并没有借助其他包。