How to limit the raster processing extent using a spatial mask?

Question

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)

Answer 1

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