用户自定义的掩膜部分输出像元的 ENVITask
12009 给本文评分:
3.0
用户自定义的掩膜部分输出像元的 ENVITask
匿名 2015年5月14日，星期四
我上一篇 Data Point 博客文章 ENVI 5.2 SP1 中的用户自定义 ENVITasks 介绍了如何使用新的 ENVITaskFromProcedure 类来创建你自己的自定义 ENVITasks，以封装任何你有的 IDL 过程。今天，我将探讨一个结合了几个更高级主题的具体示例：
- 栅格统计
- 阈值 ROI 掩膜
- 保存结果时正确掩膜像元
这个示例的灵感来源于我最近被要求协助解决的一个技术支持问题，但我对其进行了扩展，以便对读者您更有用。问题是关于如何使用 ENVIROIMaskRaster() 对栅格进行像元掩膜，然后导出结果，以便掩膜的像元可以持久保存在磁盘上。我可以只关注这部分，但我觉得展示如何使用 ENVI API 构建 ROI（而不仅仅是从文件加载）也很有用。
此任务的目标是获取一个输入栅格，评估其直方图，然后掩膜掉直方图中顶部或底部 N% 的像元。这类似于我们线性百分比拉伸模式所做的，但在那种情况下，顶部和底部的像元在其各自波段中被钳制为 0 或 255 值，中间的值则在该范围内线性缩放。此任务将使这些极端像元无效，因此它们在可视化时是透明的，并被我们的分析所忽略。
您可能担心需要遍历栅格中的所有像元来构建直方图，但别担心，ENVIRasterStatistics() 函数是在 ENVI 5.2 中引入的，它提供了直方图计算的优化实现。因此，我们可以使用此函数配合 HISTOGRAMS 和 HISTOGRAM_BINSIZE 关键字来获取每个波段的直方图。默认情况下，此函数将返回一个 256 箱的直方图，但在此任务中，我们将箱大小设置为 1，以便为每个波段获得（最大值 - 最小值 + 1）个箱。
获得直方图后，我们将其传递给带有 /CUMULATIVE 关键字的 Total()。这将产生一个与直方图大小相同的数组，其中每个元素是所有索引小于或等于输出索引的输入元素之和。每个输出数组的最后一个元素将等于栅格中的像元数。因此，如果我们除以最后一个元素（或像元数），我们将得到每个波段像元的累积分布函数 (CDF)。
然后，我们可以使用 WHERE 测试来查找有多少个箱的 CDF 值小于 N%，或者大于 (1 – N%)。这将告诉我们，如果我们想要消除顶部和底部 N% 的像元，像素值的范围将是多少。我们使用这些值来创建一个阈值 ROI，使用 ENVIROI::AddThreshold 过程。
对于如何使用 ROI 来掩膜栅格，我们有两个选择。我们可以创建 1 个 ROI 并为其每个波段添加单独的阈值，或者我们可以创建多个 ROI，每个 ROI 具有一个波段阈值。当我们使用 ENVIROIMaskRaster() 虚拟栅格函数时，结果将是相同的——只有位于 每个 阈值之外的像元才会被掩膜。这是一个逻辑与 (AND) 操作，但我想要一个逻辑或 (OR) 操作，即位于 任何 阈值之外的每个像元都将被掩膜。为此，我需要为每个波段阈值创建单独的 ROI，并使用 ENVIROIMaskRaster() 函数连续应用它们。虚拟栅格函数的一个很好的优点是它们不会将任何内容提交到磁盘，它们构建了一个栅格操作的过滤器链，这些操作基于每个瓦片对像元执行，因此这样做的性能开销非常小。
获得掩膜栅格后，我们希望将其保存到磁盘，以便掩膜的像元保持掩膜状态。我们不希望通过写出布尔值数组来浪费磁盘空间，因此使用“数据忽略值”将掩膜的像元标记为无效。不幸的是，在 ENVI 5.2 SP1 中没有一行代码的解决方案，因此我们必须为掩膜的像元“烧入”数据忽略值。在 ENVI 5.3 中，我们通过 ENVIRaster::Export 上的 DATA_IGNORE_VALUE 关键字添加了此功能，但与此同时，我们需要创建一个新的 ENVIRaster 对象并显式更改掩膜的像素值，以便我们可以持久化掩膜。因此，我们使用 ENVIRaster() 工厂函数，并使用 INHERITS_FROM 关键字从我们的输入栅格复制所有元数据和空间参考对象。我们还将提供此输出栅格将写入的 URI，并指定要使用的 DATA_IGNORE_VALUE，以便 ENVI 头文件包含该元数据字段并为我们正确掩膜像元。
ENVIRasterIterator 类使像元处理变得简单得多。我们首先需要获取一个栅格迭代器，这是通过 ENVIRaster::CreateTileIterator() 函数完成的。我们将省略关键字，因为我们希望处理整个栅格，并让它决定使用的最佳瓦片大小。使用迭代器有几种方法，虽然我通常偏爱 foreach 循环而不是 for 循环，但这次我们不能使用 foreach 语法，因为我们需要从 ENVIRasterIterator::Next() 函数获取 PIXEL_STATE 输出。一旦我们有了 PIXEL_STATE 数组，我们就使用 WHERE 查找所有非零值（已经具有数据忽略值或被掩膜的像元），然后用数据忽略值替换所有相应的像元。然后我们使用 ENVIRaster::SetTile 过程将像元复制到输出栅格，该过程使用当前的栅格迭代器对象，以便它知道像素数组的行/样本/波段范围，并将其正确提交到磁盘。一旦我们完成遍历栅格并将掩膜的像元复制到输出栅格，我们只需在输出栅格上调用 ENVIRaster::Save 即可完成将像元提交到磁盘并为我们生成头文件。
以下是我的过程的最终版本，它使用关键字来指定输入栅格、阈值百分比、数据忽略值和要使用的输出文件名：
pro maskExtremePixels, INPUT_RASTER=inputRaster, THRESHOLD_PERCENT=percent, $
DATA_IGNORE_VALUE=dataIgnoreValue, OUTPUT_FILENAME=outFilename
compile_opt idl2
; calculate raster statistics and histograms
stats = ENVIRasterStatistics(inputRaster, /HISTOGRAMS, HISTOGRAM_BINSIZE=1)
; create local variable to represent current virtual raster, starting with input raster
currRaster = inputRaster
; iterate over list of histogram hashes
foreach hist, stats['HISTOGRAMS'], band do begin
; calculate cumulative distribution function from histogram
cdf = Total(hist['COUNTS'], /Cumulative)
cdf /= cdf[-1]
; find CDF bins that are above/below threshold
lowInd = Where(cdf lt percent, lowCount)
highInd = Where(cdf gt (1.0 - percent), highCount)
; calculate digital number that corresponds to threshold bins
minVal = stats['MIN', band] + lowCount
maxVal = stats['MAX', band] - highCount
; create ROI to mask the input raster
roi = ENVIRoi(NAME='Band Thresholds', COLOR='Blue')
; add threshold definition for current band to ROI
roi.AddThreshold, inputRaster, band, MIN_VALUE=minVal, MAX_VALUE=maxVal
; apply ROI to raster to mask out extreme pixels, and update current reference
currRaster = ENVIRoiMaskRaster(currRaster, roi)
endforeach
; create raster to commit masked pixels to disk
outRaster = ENVIRaster(INHERITS_FROM=inputRaster, $
DATA_IGNORE_VALUE=dataIgnoreValue, $
URI=outFilename)
; create a raster iterator to access data
iterator = currRaster.CreateTileIterator()
for i = 1, iterator.nTiles do begin
; get the current tile's pixel and pixel value state arrays
tile = iterator.Next(PIXEL_STATE=pixelState)
; identify any pixels which should be masked out
maskInd = where(pixelState ne 0, maskCount)
; apply pixel mask if any pixels in this tile are masked
if (maskCount gt 0) then begin
tile[maskInd] = dataIgnoreValue
endif
; push masked pixels into output raster
outRaster.SetTile, tile, iterator
endfor
; save output raster and header to disk
outRaster.Save
end
现在我们想把它变成一个 ENVITask，以便我们可以在桌面和云中使用 ENVI Services Engine 运行它。我们需要做的就是构建一个描述此过程输入参数的任务模板文件：
{
"name": "MaskExtremePixels",
"baseClass": "ENVITaskFromProcedure",
"routine": "MaskExtremePixels",
"displayName": "掩膜极端像元",
"description": "此任务处理输入 ENVIRaster 上每个波段的直方图，并构建像元掩膜，排除直方图中顶部或底部百分之几的所有像元。",
"version": "5.2.1",
"parameters": [
{
"name": "INPUT_RASTER",
"displayName": "输入栅格",
"dataType": "ENVIRASTER",
"direction": "input",
"parameterType": "required",
"description": "指定要掩膜的栅格。"
},
{
"name": "THRESHOLD_PERCENT",
"displayName": "阈值百分比",
"dataType": "FLOAT",
"direction": "input",
"parameterType": "required",
"description": "从顶部和底部掩膜的直方图百分比。"
},
{
"name": "DATA_IGNORE_VALUE",
"displayName": "数据忽略值",
"dataType": "DOUBLE",
"direction": "input",
"parameterType": "required",
"description": "用于 OUTPUT_RASTER 中掩膜像元的值。"
},
{
"name": "OUTPUT_RASTER",
"keyword": "OUTPUT_FILENAME",
"displayName": "输出栅格",
"dataType": "ENVIRASTER",
"direction": "output",
"parameterType": "required",
"description": "INPUT_RASTER 的掩膜版本。"
}
]
}
在大多数情况下，这个任务模板非常直接。我们提供一个名称和描述，然后由于我们封装了一个过程，我们使用 ENVITaskFromProcedure 类作为 baseClass，并使用过程名称作为 routine。然后我们有了参数列表，包括它们的 name、displayName、dataType、direction、parameterType 和 description 属性。唯一的异常是 OUTPUT_RASTER 参数。如果我们查看过程签名，你会看到关键字 OUTPUT_FILENAME，它是一个输入字符串，告诉过程将栅格写入磁盘的位置。当我们将其包装为任务时，我们希望让任务使用者更容易，因此我们将输出参数的类型设为 ENVIRaster，并将其 keyword 设置为 OUTPUT_FILENAME。加载此任务模板时，任务框架将自动添加一个类型为 ENVIURI 的输入参数来包装此相同关键字，并将“_URI”后缀附加到原始参数的名称。因此，用户可以显式设置 OUTPUT_RASTER_URI 参数值来控制掩膜栅格的写入位置，或者他们可以不设置它，在这种情况下，任务框架将为该关键字生成一个临时文件名。
一旦你将此任务模板复制到扩展名为 .task 的文件，并将过程 PRO 文件复制到自定义代码文件夹（Windows 上为 C:\Program Files\Exelis\ENVI52\custom_code）或其他已记录的位置，你就可以使用它了。以下是一个示例，展示如何使用：
e = envi()
file = FilePath('qb_boulder_msi', ROOT_DIR=e.ROOT_DIR, SUBDIR='data')
oRaster = e.OpenRaster(file)
oTask = enviTask('maskExtremePixels')
oTask.INPUT_RASTER = oRaster
oTask.THRESHOLD_PERCENT = 0.02
oTask.DATA_IGNORE_VALUE = -1
oTask.Execute
oView = e.getView()
oLayer = oView.CreateLayer(oTask.OUTPUT_RASTER)