如何使用空间掩码限制栅格处理范围？

Question

我正在尝试将 MATLAB 中的栅格处理限制为仅包含 shapefile 边界内的区域，类似于 ArcGIS Spatial Analyst 函数如何使用 mask。这是我正在使用的一些（可重现的）示例数据：

• 4-band NAIP image（警告 169MB 下载）
• shapefile of study area boundaries（File Dropper 上的压缩 shapefile）

这是我用来计算 NDVI 的 MATLAB 脚本：

file = 'C:\path\to\doi1m2011_41111h4nw_usda.tif';
[I R] = geotiffread(file);
outputdir = 'C:\output\'

% Calculate NDVI
NIR = im2single(I(:,:,4));
red = im2single(I(:,:,1));

ndvi = (NIR - red) ./ (NIR + red);
double(ndvi);
imshow(ndvi,'DisplayRange',[-1 1]);

% Stretch to 0 - 255 and convert to 8-bit unsigned integer
ndvi = floor((ndvi + 1) * 128); % [-1 1] -> [0 256]
ndvi(ndvi < 0) = 0;             % not really necessary, just in case & for symmetry
ndvi(ndvi > 255) = 255;         % in case the original value was exactly 1
ndvi = uint8(ndvi);             % change data type from double to uint8

% Write NDVI to .tif file (optional)
tiffdata = geotiffinfo(file);
outfilename = [outputdir 'ndvi_' 'temp' '.tif'];  
geotiffwrite(outfilename, ndvi, R, 'GeoKeyDirectoryTag', tiffdata.GeoTIFFTags.GeoKeyDirectoryTag) 

下图说明了我想使用 MATLAB 完成的工作。对于此示例，我使用 ArcGIS raster calculator (Float(Band4-Band1)/Float(Band4+Band1)) 来生成右侧的 NDVI。我还将研究区域 shapefile 指定为 mask in the environment settings。

问题：

如何在 MATLAB 中使用多边形 shapefile 作为空间掩码来限制栅格处理范围以复制图中所示的结果？

我不成功尝试了什么：

roipoly 和 poly2mask，尽管我似乎无法正确应用这些功能（考虑到这些是空间数据）来产生所需的效果。

编辑：

我尝试了以下方法将 shapefile 转换为蒙版，但没有成功。不知道我哪里出错了......

s = 'C:\path\to\studyArea.shp'

shp = shaperead(s)
lat = [shp.X];
lon = [shp.Y];

x = shp.BoundingBox(2) - shp.BoundingBox(1)
y = shp.BoundingBox(3) - shp.BoundingBox(1) 

x = poly2mask(lat,lon, x, y)

---

错误信息：

Error using poly2mask
Expected input number 1, X, to be finite.

Error in poly2mask (line 49)
validateattributes(x,{'double'},{'real','vector','finite'},mfilename,'X',1);

Error in createMask (line 13)
x = poly2mask(lat,lon, x, y)

Answer 1

这里分为三个步骤，我将为此创建3个函数：

1. 计算完整输入图像的 NDVI：ndvi = comp_ndvi(nir, red)
2. 从 shapefile 计算掩码：mask = comp_mask(shape)
3. 结合 NDVI 和掩码：output = combine_ndvi_mask(ndvi, mask)

您的问题中有comp_ndvi() 的代码。 combine_ndvi_mask() 的代码取决于您要对蒙版区域执行的操作；如果你想让它们变成白色，它可能看起来像：

function output = combine_ndvi_mask(ndvi, mask)
output = ndvi;
output(~mask) = 255;
end

在comp_mask() 中，您需要使用poly2mask() 将多边形顶点转换为光栅蒙版。为了在这里提供帮助，我需要知道你已经拥有了什么。您是否已将顶点加载到 MATLAB 中？你用 poly2mask 试过什么？

Answer 2

您可以通过以下方式读取感兴趣的区域：

roi = shaperead('study_area_shapefile/studyArea.shp');

砍掉尾随的 NaN：

rx = roi.X(1:end-1);
ry = roi.Y(1:end-1);

如果您的 shapefile 中有多个多边形，它们由 NaN 分隔，您必须单独处理它们。

然后使用来自卫星图像空间参考的 worldToIntrinsic 方法将多边形点转换为图像坐标：

[ix, iy] = R.worldToIntrinsic(rx,ry);

这假设两个坐标系相同。

然后你可以去制作你的面具：

mask = poly2mask(ix,iy,R.RasterSize(1),R.RasterSize(2));

在进行任何计算之前，您可以在原始多层图像上使用蒙版：

I(repmat(~mask,[1,1,4])) = nan;

或通过以下方式在单层（即红色）上使用它：

red(~mask) = nan;

如果区域非常小，将蒙版图像转换为稀疏矩阵可能是有益的（对于内存和计算能力）。我还没有尝试过这是否会产生任何速度差异。

red(~mask) = 0;
sred = sparse(double(red));

不幸的是，稀疏矩阵只能使用双精度数，因此您的 uint8 需要在转换之前。

通常，您应该从图像中裁剪出 ROI。查看对象“roi”和“R”以找到有用的参数和方法。我这里没做过。

最后是我的脚本版本，稍有其他改动：

file = 'doi1m2011_41111h4nw_usda.tif';
[I R] = geotiffread(file);
outputdir = '';

% Read Region of Interest
roi = shaperead('study_area_shapefile/studyArea.shp');
% Remove trailing nan from shapefile
rx = roi.X(1:end-1);
ry = roi.Y(1:end-1);
% convert to image coordinates
[ix, iy] = R.worldToIntrinsic(rx,ry);
% make the mask
mask = poly2mask(ix,iy,R.RasterSize(1),R.RasterSize(2));
% mask sat-image
I(repmat(~mask,[1,1,4])) = 0;

% convert to sparse matrizes
NIR = sparse(double(I(:,:,4)));
red = sparse(double(I(:,:,1)));
% Calculate NDVI
ndvi = (NIR - red) ./ (NIR + red);
% convert back to full matrizes
ndvi = full(ndvi);
imshow(ndvi,'DisplayRange',[-1 1]);

% Stretch to 0 - 255 and convert to 8-bit unsigned integer
ndvi = (ndvi + 1) / 2 * 255; % [-1 1] -> [0 255]
ndvi = uint8(ndvi);          % change and round data type from double to uint8 

% Write NDVI to .tif file (optional)
tiffdata = geotiffinfo(file);
outfilename = [outputdir 'ndvi_' 'temp' '.tif'];  
geotiffwrite(outfilename, ndvi, R, 'GeoKeyDirectoryTag', tiffdata.GeoTIFFTags.GeoKeyDirectoryTag);
mapshow(outfilename);