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

Answer 2

There are three steps here, for which I will create 3 functions:

1. Compute the NDVI for the complete input image: ndvi = comp_ndvi(nir, red)
2. Compute the mask from the shapefile: mask = comp_mask(shape)
3. Combine the NDVI and the mask: output = combine_ndvi_mask(ndvi, mask)

You have the code for comp_ndvi() in your question. The code for combine_ndvi_mask() depends on what you want to do to the masked areas; if you want to make them white, it might look like:

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

In comp_mask() you will want to use poly2mask() to convert the polygon vertices into the raster mask. In order to help here I need to know what you've got already. Have you loaded the vertices into MATLAB? What have you tried with poly2mask?