Как ограничить размер растровой обработки с помощью пространственной маски?

Question

Я пытаюсь ограничить растровую обработку в MATLAB, чтобы включить только области внутри границы шейп-файла, аналогично тому, как функции ArcGIS Spatial Analyst используют mask. Вот некоторые (воспроизводимые) выборка данных Я работаю с:

• 4-band NAIP image (ВНИМАНИЕ 169MB скачать)
• shapefile of study area boundaries (заархивированный шейпфайл на File Dropper)

Вот скрипт MATLAB Я использую для расчета 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) 

на следующем рисунке показано, что я хотел бы достичь с помощью MATLAB. В этом примере я использовал ArcGIS raster calculator (Float (Band4-Band1)/Float (Band4 + Band1)) для создания NDVI справа. Я также указал форма формы области исследования как mask in the environment settings.

Вопрос:

Как можно ограничить степень растровой обработки в MATLAB, используя многоугольник шейп в качестве пространственной маски, чтобы повторить результаты, показанные на рисунке?

То, что я безуспешно пытался:

roipoly и poly2mask, хотя я не могу применить эти функции должным образом (с учетом их пространственные данные) для получения желаемого эффекта.

EDIT:

Я попытался следующие преобразовать шейп в маску, но безуспешно. Не знаете, где я неправильно здесь ...

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

Вы можете прочитать интересующую область по:

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

Чоп трейлинг NaN:

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

Если у вас есть несколько многоугольников в вашем шейп, они разделяются NaNs и у вас есть рассматривать их отдельно.

Затем используют 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)); 

Unfortunatly, разреженные matrizes возможны только с двойниками, так что ваш 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); 

Answer 2

Есть три шага здесь, для которого я буду создавать 3 функции:

1. Вычислить NDVI для полного входного изображения: ndvi = comp_ndvi(nir, red)
2. Вычислить m просить у шейп: 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?