如何使用空间掩码限制栅格处理范围？ [英] How to limit the raster processing extent using a spatial mask?
问题描述
我正在尝试限制MATLAB中的栅格处理以仅包含shapefile边界内的区域，类似于ArcGIS Spatial Analyst函数如何使用和 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）
错误消息：
使用poly2mask时出错
预期输入数1，X为有限的。
poly2mask中的错误（第49行）
validateattributes（x，{'double'}，{'real'，'vector'，'finite'}，mfilename，'X'，1 ）;
createMask中的错误（第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分开，你必须单独对待它们。
然后使用sat-image空间参考中的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;
如果区域非常小，转换a可能是有益的（对于内存和计算能力）掩蔽图像到稀疏矩阵。我没试过，如果它有任何速度差异。
red（〜mask）= 0 ;
sred =稀疏（双（红色））;
不幸的是，稀疏的matrizes只能用双打，所以你的uint8需要在转换之前。
通常你应该从图像中裁剪出ROI。查看对象roi和R以查找有用的参数和方法。我还没有在这里完成。
最后我的脚本版本，稍作其他一些改动：
file ='doi1m2011_41111h4nw_usda.tif';
[I R] = geotiffread（file）;
outputdir ='';
％阅读感兴趣的区域
roi = shaperead（'study_area_shapefile / studyArea.shp'）;
％从shapefile中删除尾随nan
rx = roi.X（1：end-1）;
ry = roi.Y（1：end-1）;
％转换为图像坐标
[ix，iy] = R.worldToIntrinsic（rx，ry）;
％使掩码
掩码= poly2mask（ix，iy，R.RasterSize（1），R.RasterSize（2））;
％mask sat-image
I（repmat（~mask，[1,1,4]））= 0;
％转换为稀疏matrizes
NIR =稀疏（double（I（：，：，4）））;
red =稀疏（double（I（：，：，1）））;
％计算NDVI
ndvi =（NIR - 红色）./（NIR +红色）;
％转换回完整matrizes
ndvi = full（ndvi）;
imshow（ndvi，'DisplayRange'，[ - 1 1]）;
％拉伸至0 - 255并转换为8位无符号整数
ndvi =（ndvi + 1）/ 2 * 255; ％[-1 1] - > [0 255]
ndvi = uint8（ndvi）; ％更改和舍入数据类型从double到uint8
％将NDVI写入.tif文件（可选）
tiffdata = geotiffinfo（file）;
outfilename = [outputdir'ndvi_''temp''.tif'];
geotiffwrite（outfilename，ndvi，R，'GeoKeyDirectoryTag'，tiffdata.GeoTIFFTags.GeoKeyDirectoryTag）;
mapshow（outfilename）;
I am trying to limit raster processing in MATLAB to include only areas within a shapefile boundary, similar to how ArcGIS Spatial Analyst functions use a mask. Here is some (reproducible) sample data I am working with:
- A 4-band NAIP image (WARNING 169MB download)
- A shapefile of study area boundaries (A zipped shapefile on File Dropper)
Here is a MATLAB script I use to calculate NDVI:
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)
The following image illustrates what I would like to accomplish using MATLAB. For this example, I used the ArcGIS raster calculator (Float(Band4-Band1)/Float(Band4+Band1)) to produce the NDVI on the right. I also specified the study area shapefile as a mask in the environment settings.
Question:
How can I limit the raster processing extent in MATLAB using a polygon shapefile as a spatial mask to replicate the results shown in the figure?
What I have unsuccessfully tried:
roipoly and poly2mask, although I cannot seem to apply these functions properly (taking into account these are spatial data) to produce the desired effects.
EDIT:
I tried the following to convert the shapefile to a mask, without success. Not sure where I am going wrong here...
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 messages:
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)
You can read the region of interest by:
roi = shaperead('study_area_shapefile/studyArea.shp');
Chop the trailing NaN:
rx = roi.X(1:end-1);
ry = roi.Y(1:end-1);
If you have several polygons in your shapefile, they are seperated by NaNs and you have to treat them seperately.
Then use the worldToIntrinsic-method from your spatial reference of the sat-image to convert the polygon-points into image-coordinates:
[ix, iy] = R.worldToIntrinsic(rx,ry);
This assumes both coordinate systems are the same.
Then you can go and make your mask by:
mask = poly2mask(ix,iy,R.RasterSize(1),R.RasterSize(2));
You can use the mask on your original multilayer image before making any calculation by:
I(repmat(~mask,[1,1,4])) = nan;
Or use it on a single layer (i.e. red) by:
red(~mask) = nan;
If the regions are very small, it could be beneficial (for memory and computation power) to convert a masked image to a sparse matrix. I have not tried if that makes any speed-difference.
red(~mask) = 0;
sred = sparse(double(red));
Unfortunatly, sparse matrizes are only possible with doubles, so your uint8 needs prior to be converted.
Generally you should crop the ROI out of the image. Look in the objects "roi" and "R" to find useful parameters and methods. I haven't done it here.
Finally my version of your script, with some slight other changes:
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);
这篇关于如何使用空间掩码限制栅格处理范围？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
