基于Python的Landsat影像建设用地自动识别与提取

 

 

 

    系统实现基于Python语言进行编码，数据预处理和质量控制部分通过自行编码实现，主要空间分析函数调用ArcGIS中的空间分析模块。

 

1 数据预处理

 

    (1)对下载的原始数据进行批量解压。对海量原始影像压缩文件夹进行遍历，识别并逐次解压。文件解压关键代码如下:

 

targzfilc. cxtract(filc,  unzipcdfoldcr)#从压缩文件targzfile中把解压出的文件[file放到文件夹nzipcdfoldcr里。

 

    (2)重命名。在系统中，为了后续流程的完整性，对已解压文件名不统一的波段重新命名。文件重命名关键代码如下:

    arcpy. Rcnamc_managcmcnt(namc_fr,  name _to)#从将文件name_fr重命名为name_to。

    (3)提取关键波段。海量影像的所有波段中，依据建设用地识别方法，对波段进行筛选，单独提取出关键的波段数，并进行另存。

    (4)获取最小重叠区。去除空值和影像边界，其留下的最小边界作为后面研究范围。以关键波段中某一波段为代表，利用Python脚本对其影像进行批量化掩摸，计算最小重叠区。关键代码如下:

    arcpy. gp.CelllStatistics  sa_(in,out,”MINIMUM";”DATA" )#获取所有输人图层的最小值。剔除外围以值表不的空值黑边，只保留有值区，裁剪出最小重叠区。关键代码如下:

    arcpy, gp. Con_sa(in，1 ， out,”，”Value>0" )#剔除外围黑边。

 

2 云噪声去除

 

为了消除时相干扰，去除云噪声，需进行时域中值滤波。其原理是云的反射率在各波段都位于异常极大值区，而云影大多位于异常极小值区，取中值具有很好的去云和去云影效果。关键代码如下：

    arcpy. gp. CellStatistics _sa(m,  out,”MEDIAN"，”DATA'. )#获取所有输人图层的中值。

 

3 稳定像元获取

 

    不同影像中波段值存在时相差异，需对生成的波段中值进行z值标准化，并计算影像的z值标准差，再以0.5倍标准差为闽值，替换不稳定像元，最后反z值标准化，即基于统一的平均值和标准差系数反算回去，获取每个单波段的稳定像元图层。该步骤能保证在后续步骤中影像均采用相同阈值进行分割。关键代码如下:

    arcpy. gp. Minus_sa(in, MEAN, numerator)#获得每期影像与多期影像平均值之差，以作为标准化的分子numerator

    arcpy. gp. Divide_sa(numerator, STD, out)#获得每期影像与numerator之比，得到该期影像的标准化值

    arcpy, gp. Con_sa(in , 1 , out,”，”Value <0. 5" )#获得均值两侧0.5倍标准差范围内的像元，即稳定像元。

 

    (1)根据波段计算出VDVI,MVDWI指数。

    (2)采用归一化闽值法确定植被的闽值范围为0.25 ,粗略去除植被。

    arcpy. gp. Con_sa<(MNDWI，1，Pout，0，'Value>0. 25')#获取像元大于0. 25的NDVI，用于识别植被。

    C3)采用步骤C2)中的归一化闽值法，确定水体的闽值范围为0. 1 ,去除水体。

    arcpy. gp. Con_sa<(MNDWI，1，Pout，0，'Value>0.10')#获取像元大于0.1的MNDWI，用于识别水体

    (4)在前3步的基础上，通过光谱和纹理分析，去除试验区内比较明显的植被、水体以后，通过光谱差异和借用实验区范围内的Google Earth影像对比，进一步识别出需要识别去除的高亮地物(含裸土、部分滩涂、部分高亮人工地面、部分高亮屋顶)、水体及偏阴坡植被、裸地及植被。

 

    arcpy, gp. Con_sa(D73，1，Pout，0，'Value>0. 06')#获取像元大于0. 06的D73，用于识别高亮地物

    arcpy. gp. Con_sa(D73 , 1，Pout,O,'Valuc<0. O 1')#获取像元小于0. 0l的D73，用于识别水体及偏阴坡植被

    arcpy. gp. Con_sa ( ND75ND42，1，Pout，0，' Value<一0.40')#获取像元大于-0. 04的ND75ND42，识别裸地及植被。

 

4 建设用地精化提取

 

    初步识别可能的建设用地后，系统用于精化提取建设用地的方法主要是采用空间邻域关系计算。分别对水体、植被的边界像元通过邻域均值计算后，采用条件函数逐一判别消除其它地物。其关键代码如下:

 

    arcpy, gp. Con_sa<(Vm，1，Vm01，0，”VALUE>0. 25")#获取邻域内均值Vm超过0. 25的植被

    arcpy, gp. Con_sa<(Wm，1，W m01，0，”VALUE>0.25" )##获取邻域内均值Wm超过0. 25的水体

    arcpy.gp.CellStatistics_sa(Vm01&Wm01，Vm01 sumWm01 , " SUM" , " DATA" )#获取邻域内水体和植被Vm01g&Wm01均超过0. 25的像元的和

    arcpy, gp. Con_sa<(VmOlsumWm01,1,Pout，0，”VALUE=2”)#获取上一步中VALUE=2的值。

 

 

                               建设用地提取结果

 

更多遥感知识，请关注我的微信公众号，扫下二维码即可关注。