将分类像素转换为 Shapefile
17234 给本文评分:
5.0
将分类像素转换为 Shapefile
匿名作者 2016年5月19日，星期四
通常，当我在处理一个项目时，最终会生成一幅图像作为结果，这很好。分析人员可以查看图像，并据此做出决策。然而，也有很多时候，我的结果需要输入到另一个程序中进行进一步处理。这种交接可以通过多种不同的方式完成，但最常见的方式之一是使用 shapefile。唯一的问题是，如何将所需的数据从栅格空间转换到矢量空间？
下面，我设想了一个场景，需要提取湖泊的边界。我的初始输入是一幅 Landsat 8 影像，以及一个使用 ROI 工具绘制的感兴趣区域（ROI）。我的目标是生成一个简单的 shapefile，其中包含下图所示影像中的湖岸线。在处理单个像素时，这并不是一项非常困难的任务，但在将像素转换为矢量时，可能会变得有些棘手。
图 1：带有绿色 ROI 叠加层的原始影像
为了开始这个过程，我将所有数据导入 IDL，并获取当前 ENVI 实例的引用。
; 获取当前 envi 实例 e = ENVI(/current)
; 获取数据集合 DataColl = e.Data
; Landsat 8 OLI 数据集每个波段对应一个 TIFF 文件， ; 并有一个关联的元数据文件 (*_MTL.txt)。打开 ; 元数据文件可以自动读取增益和偏移量。
IMG_File = File_Search('C:\MyData, '*_MTL.txt') Raster = e.OpenRaster(IMG_File)
; 一个包含研究区域的 ENVI ROI ROI_File = File_Search('C:\MyData\ROIs', '*.xml') oROI = e.OpenROI(ROI_File[0], PARENT_RASTER=Raster[0])
; 设置输出 shapefile 的名称 output_SHP = 'C:\Output\Lake.shp'
将数据导入 IDL 后，我对图像进行了裁剪，使其仅包含我的感兴趣区域，并将数字量化值转换为辐亮度值。
; 获取 ROI 的范围 extent = oROI.GetExtent(Raster[0])
; 将图像裁剪到研究区域的范围 Subset = OBJARR(n_elements(Raster)) foreach Band, Raster, index do Subset[index] = Band.Subset(Sub_rect = extent)
; 从 ENVI 任务目录中获取辐射定标任务 RadCal_Task = ENVITask('RadiometricCalibration')
; 定义输入。由于辐亮度是默认的校准方法， ; 此处无需指定。 RadCal_Task.Input_Raster = Subset[0] ; Bands 1-7 RadCal_Task.Output_Data_Type = 'Double'
; 定义输出栅格 URI RadCal_Task.Output_Raster_URI = e.GetTemporaryFilename()
; 运行任务 RadCal_Task.Execute
现在影像数据处于可用状态，我计算了归一化差异水体指数（NDWI），并对结果进行阈值分割，以提取大面积水体。
; 计算修正的归一化差异水体指数（NDWI）
; Xu, H. Modification of normalized difference water index (NDWI) to enhance ; open water features in remotely sensed imagery. Int. J. Remote Sens. 2006, ; 27, 3025–3033.
Green = RadCal_Task.Output_Raster.GetData(Bands=2) SWIR = RadCal_Task.Output_Raster.GetData(Bands=5)
MNDWI = (Green - SWIR) / (Green + SWIR)
; 制作仅包含水体的二值掩膜 Water = MNDWI gt .7
图 2：水体区域的二值掩膜
由于我只对图像中心的大湖泊感兴趣，我找到哪个区域包含最多的像素，并将其分离出来。
; 查找图像中有多少个不同的区域 Labeled_Water = LABEL_REGION(water,/ALL_NEIGHBORS)
; 查找每个区域的像素数（排除零值区域） Region_Stats = lonarr(max(Labeled_Water)) foreach region, findgen(max(Labeled_Water)-1)+1, index do $ Region_Stats[index] = n_elements(where(Labeled_Water eq region))
; 选择像素最多的区域 junk = max(Region_Stats,Region_pos)
; 将零位置加回来 Region_pos++
通过使用 CONTOUR 函数，可以找到湖泊边缘的所有像素。这在我们开始构建 shapefile 时非常重要。
; 获取目标边界周围所有像素的位置 Labeled_Water = Labeled_Water eq Region_pos contour, Labeled_Water, PATH_INFO=path_info, path_xy=path_xy,$ /PATH_DATA_COORDS, LEVELS=[1]
; 确保第一个点和最后一个点匹配 if path_xy[0,-1] ne path_xy[0,0] and path_xy[1,-1] ne path_xy[1,0] then $ path_xy = [[path_xy], [path_xy[*,0]]]
; 从文件坐标转换到地图坐标 SpatialRef = RadCal_Task.Output_Raster.SpatialRef SpatialRef.ConvertFileToMap, path_xy[0,*],path_xy[1,*], MapX, MapY
现在所有湖岸线像素都已处理完毕，可以构建最终的输出——shapefile 了。
; 创建 shp 文件 oShp = obj_new('IDLffShape', output_SHP, /UPDATE, ENTITY_TYPE=5)
; 添加属性 oShp->AddAttribute, 'Name', 7, 15
; 创建新的实体 ent = {IDL_SHAPE_ENTITY}
; 设置实体索引 entity_index = 0
; 填充实体信息 ent.SHAPE_TYPE = 5 ent.ISHAPE = entity_index ent.N_VERTICES = n_elements(mapX) ent.BOUNDS[0] = min(mapx) ent.BOUNDS[1] = min(mapy) ent.BOUNDS[2] = 0.0 ent.BOUNDS[3] = 0.0 ent.BOUNDS[4] = max(mapx) ent.BOUNDS[5] = max(mapy) ent.BOUNDS[6] = 0.0 ent.BOUNDS[7] = 0.0 ent.VERTICES = ptr_new([MapX,MapY])
; 填充属性 attr = oShp->GetAttributes(/ATTRIBUTE_STRUCTURE) attr.ATTRIBUTE_0 = 'My Target'
; 将实体添加到 shapefile oShp->PutEntity, ent
; 将属性添加到 shapefile oShp->SetAttributes, entity_index, attr
; 关闭 shapefile OBJ_DESTROY, oShp
; 创建投影文件 prj_name = file_dirname(output_SHP) + path_sep() + file_basename(output_SHP, '.shp') + '.prj' openw, lun, prj_name, /get_lun writeu, lun, byte(Spatialref.COORD_SYS_STR) close, lun free_lun, lun
图 3：显示最终 shapefile（红色）的影像
我希望这有助于解释如何将几个不同的 IDL 和 ENVI 功能链接起来，以实现从栅格到矢量的转换。