如何使用空间掩模限制光栅处理范围？

我试图将 MATLAB 中的栅格处理限制为仅包含 shapefile 边界内的区域，类似于 ArcGIS Spatial Analyst 函数使用mask http://resources.arcgis.com/en/help/main/10.2/index.html#//001w0000001t000000。这是我正在使用的一些（可重现的）示例数据：

• A 4 波段 NAIP 图像 http://cloud.insideidaho.org/data/imageryBaseMapsEarthCover/orthoimagery/2011_1m_idaho/41/doi1m2011_41111h4nw_usda.tif（警告下载 169MB）
• A 研究区域边界的shape文件 http://www.filedropper.com/studyareashapefile（File Dropper 上的压缩 shapefile）

这是我用来计算的 MATLAB 脚本NDVI http://resources.arcgis.com/en/help/main/10.2/index.html#/NDVI_function/009t00000052000000/:

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 栅格计算器 http://resources.arcgis.com/en/help/main/10.2/index.html#//009z000000z7000000（Float(Band 4-Band 1)/Float(Band 4+Band 1)) 生成右侧的 NDVI。我还将研究区域形状文件指定为环境设置中的掩码 http://resources.arcgis.com/en/help/main/10.2/index.html#//001w0000001t000000.

问题：

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

我拥有的不成功 tried:

roipoly http://www.mathworks.com/help/images/ref/roipoly.html and 多聚二掩模 http://www.mathworks.com/help/images/ref/poly2mask.html，尽管我似乎无法正确应用这些函数（考虑到这些是空间数据）来产生所需的效果。

EDIT:

我尝试了以下方法将 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)

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

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);