python提取二值栅格上边界和中线

碎碎念：
比较简单粗暴的方法…缺点就是提取上边界的时候有一些细碎的像元还是没办法删去，如果图案过分扭曲，效果就差很多。（初学python…代码比较冗长繁琐[扶额]）
适用于值为0和1的栅格图像。

1.提取上边界

原图像：

提取后：

代码：

def Get_UpperBoundary(dataarray,rows,cols):
    boundary=dataarray.copy()
    boundary[boundary > 0] = 0
    for i in range(cols-1):
        for j in range(rows-1):
            if (dataarray[j, i] != 0) and (((dataarray[j-1, i] == 0) and (dataarray[j,i-1]==0)) or ((dataarray[j-1,i]==0) and (dataarray[j,i+1]==0)) or ((dataarray[j-1,i]==0) and (dataarray[j-1,i-1]==0))or ((dataarray[j-1,i]==0) and (dataarray[j-1,i+1]==0))):
                boundary[j, i] = 1

    for i in range(cols-1):
        for j in range(rows-1):
           if (dataarray[j, i] != 0):
                boundary[j, i] = 1
                break
    #删去底部多余的像元点
    for i in range(cols-1):
        a = 0; b = []; c = []
        c = np.array(c)
        # 提取出非0点的位置
        for j in range(rows-1):
            if(boundary[j, i] != 0):
                b.append(j)
                a = a+1
        if(a>1):
            for k in range(len(b) - 1):
                a = b[k + 1] - b[k]  # 后者减前者
                c = np.append(c, a)  # 添加元素到新列表
            #差值大于5像元之后的其他像元会被剔除
            c[c > 5] = 999;t=0
            for n in range(len(c)):
                if(c[n] == 999):
                    t += 1;a = int(n+1)
                    break
            if(t!=0):
                aa = int(len(b)-1)
                if(a == aa):
                    boundary[b[a], i] = 0
                for m in range(a, aa):
                    boundary[b[m], i] = 0
    return boundary

2.提取中线

原图像同上
提取后：

代码：

def Get_MedianBoundary(dataarray,rows,cols):
    boundary=dataarray.copy()
    boundary[boundary > 0] = 0
    for i in range(cols-1):
        b = 0
        d = 0
        for j in range(rows-1):
            if (dataarray[j,i]!=0) :
                b += j
                d = d + 1
        if(d!=0):
            boundary[b//d, i] = 1
    return boundary

---

11.18 更新（转为连续的矢量线）

就关于线不连续的问题，可以使用arcpy进行解决，但是前提是得安装了arcgis，最终的结果会是一条连续的矢量线。使用arcpy是为了方便批量处理，如果只有一两张待处理数据，直接用arcgis做完全可以。

主要思路如下：
1.首先将栅格转点
2.选出点值为1的点要素（结果如下图所示，图中使用的为上面提取的中线数据）

3.添加字段，计算这些点的经度值（做这一步是因为arcgis在将点连成线的时候，默认是一行一行的连的，结果就会变成下图这样）
这是错误示例：

4.添加字段后，根据计算出的经度的顺序，连接点成矢量线

接下来贴上arcpy代码：

# -*- coding: utf-8 -*-
import os
import arcpy

arcpy.env.workspace =""#这里是待处理数据的文件夹位置
path = ""#待处理数据位置
PointPath = ""#栅格转点的暂存位置，后面会删除
SelectPointPath = ""#筛选点的暂存位置，后面会删除
OutpolylinePath = ""#最终结果的位置
arcpy.RasterToPoint_conversion (path, PointPath)#栅格转点
arcpy.Select_analysis(PointPath, SelectPointPath,'"GRID_CODE" = 1')#筛选点
arcpy.AddField_management(SelectPointPath, "COLS", "DOUBLE", "", "", "#", "#", "#", "#")#添加字段
arcpy.CalculateField_management(SelectPointPath, "COLS", "!shape.centroid.x!", "PYTHON_9.3")#计算各点经度
arcpy.PointsToLine_management (SelectPointPath, OutpolylinePath, "GRID_CODE", "COLS", "")#点连成线
os.remove(PointPath)
os.remove(SelectPointPath)#删除暂存的点文件