Crop and Mask Large GeoTIFF File Using Shapefile

Since R2023b

This example shows how to crop, mask, and export a large GeoTIFF file using a shape from a shapefile.

The goals of this example are to crop a large GeoTIFF file to the smallest extent that encompasses the shape in the shapefile, and to highlight the image data within the shape by applying a mask that sets the value of all pixels outside the shape to 0. Applying the cropping and masking operations directly to a large file requires a lot of memory and potentially wastes computational resources. The example mitigates these challenges by working with the GeoTIFF file as a blocked image. The blocked image approach can scale to very large files because it loads and processes one block of data at a time. The blocked image approach is also computationally efficient because it can omit blocks of nonmeaningful data, such as regions outside the shape.

The steps to process a large GeoTIFF file with cropping and masking operations include:

• Read the GeoTIFF file into a blockedImage object.

• Load the shape from the shapefile.

• Calculate the bounding box of the shape in raster coordinates, and crop the high-resolution image to the bounding box.

• Convert the shape from a set of polygon coordinates to a low-resolution raster mask.

• Mask blocks of the cropped high-resolution image that overlap with the low-resolution mask of the shape.

• Write the processed data to a TIFF file, and insert the required GeoTIFF tags to convert the file to a GeoTIFF file.

Create Multiresolution Blocked Image

Create a blockedImage object that reads the GeoTIFF file. This image is of Boston and the surrounding region. The image has one resolution level.

bim = blockedImage("boston.tif");

Create an adapter that reads TIFF files and applies JPEG compression.

adapter = images.blocked.TIFF;
adapter.Compression = Tiff.Compression.JPEG;

Create a two-level blockedImage object by using the makeMultiLevel2D function. Specify the block size 1024-by-1024 pixels, which is generally appropriate for various image sizes and machine configurations. Specify the scale factors as the two-element vector [1 1/8]. The image at the first resolution level consists of the original large image in bim. The image at the second resolution level is 1/8 the size of the original image in each dimension.

timestamp = datetime("now",Format="yyyy-MM-dd-HH-mm-ss");
outputFile = "boston_multiLevel_" + string(timestamp) + ".tif";
bim2 = makeMultiLevel2D(bim,BlockSize=[1024 1024],Scales=[1 1/8], ...
    OutputLocation=outputFile,Adapter=adapter);

Display the multiresolution blocked image in raster coordinates by using the bigimageshow function. As you zoom in and out of the displayed image, the bigimageshow function switches between the resolution levels.

bigimageshow(bim2)

GeoTIFF files have raster reference information embedded in the metadata. This information is required to transform the detection of an object in pixel coordinates to map coordinates. Get the raster reference object for the high-resolution image using the georasterinfo function.

ginfo = georasterinfo(bim.Source);
R = ginfo.RasterReference;

Load Shapefile

Read the shape from the shapefile by using the readgeotable function. The shapefile contains a hand drawn outline of two parks in the city of Boston. The shape uses the same projected coordinate reference system (CRS) as the image.

shptbl = readgeotable("bostonparks.shp");

To display the shape in map coordinates, load the low-resolution image from the blocked image. Update the corresponding raster reference object to match the size of this level.

imLowRes = gather(bim2);
RLowRes = R;
RLowRes.RasterSize = size(imLowRes);

Display the shape over the low-resolution image in map coordinates by using the mapshow function.

figure
mapshow(imLowRes,RLowRes)
mapshow(shptbl,EdgeColor="y",FaceColor="none")

Crop Image to Extents of Shape

Crop the image to the smallest extent that encompasses the shape.

Convert the coordinates of the shape from map coordinates to raster coordinates by using the raster reference object.

tbl = geotable2table(shptbl,["XWorld" "YWorld"]);
xWorld = tbl.XWorld{1};
yWorld = tbl.YWorld{1};
[x,y] = worldToIntrinsic(R,xWorld,yWorld);
parkPoly = [x(:) y(:)];

Find the bounding box of the shape in raster coordinates.

[xmin,xmax] = bounds(parkPoly(:,1));
[ymin,ymax] = bounds(parkPoly(:,2));
cropStart = floor([ymin xmin]);
cropStart = max(cropStart,1);
cropEnd = ceil([ymax xmax]);
cropEnd = min(cropEnd,bim2.Size(1,1:2));

Crop the blocked image to the bounding box using the crop function. The crop function updates the blocked image metadata with the new extents of the image at each resolution level.

bimCrop = crop(bim2,cropStart,cropEnd);

Update the corresponding reference object to use the extents of the cropped image.

rowLimits = [cropStart(1) cropEnd(1)];
colLimits = [cropStart(2) cropEnd(2)];
Rcrop = cropToBlock(R,rowLimits,colLimits);

To verify the results of the cropping operation, display the shape as a polygon over the cropped image. The bigimageshow and drawpolygon functions display the cropped blocked image and the shape polygon using raster coordinates.

bigimageshow(bimCrop)
hp = drawpolygon(Position=parkPoly);

Create Raster Mask

You can process large images more quickly by applying the operations to only regions of interest (ROIs). For this example, the ROI consists of pixels within the shape polygon. Convert the polygon vertices to a raster mask at the coarse resolution level by using the poly2BlockedImage function. Ensure that the first block is full by specifying a block size of 256-by-256 pixels.

lvl = bim2.NumLevels; 
maskSize = bim2.Size(lvl,[1 2]);
mask = polyToBlockedImage({parkPoly},true,maskSize,BlockSize=[256 256], ...
    WorldStart=bimCrop.WorldStart(1,:),WorldEnd=bimCrop.WorldEnd(1,:));

This example performs the optional step of dilating the mask. Dilation extends the masked region slightly beyond the initial shape boundary. This step can improve the result of masking the image when you know the shape boundary to be slightly inaccurate around small boundary features, or you need to include additional image details just past the initial boundary.

Dilate the mask at the coarse resolution level by using the imdilate function. Call the imdilate function on all blocks in the mask by using the apply function.

mask = apply(mask,@(block)imdilate(block.Data,ones(3)));

Display the mask in raster coordinates by using the bigimageshow function.

bigimageshow(mask)

Apply Raster Mask to Cropped Image

Identify the image blocks at the high resolution level that are within the ROI by using the selectBlockLocations function and specifying the Masks name-value argument. To include all image blocks with at least one pixel within the ROI, also specify the InclusionThreshold name-value argument as 0. Ensure that the first block is full by specifying a block size of 256-by-256 pixels.

bls = selectBlockLocations(bimCrop,BlockSize=[256 256], ...
    Masks=mask,InclusionThreshold=0);

Apply the mask to the image at the finest resolution level by using the applyMask helper function, which is defined at the end of the example. Call the applyMask helper function on all blocks of the high-resolution image by using the apply function. To limit processing to the identified blocks within the ROI, specify the BlockLocationSet name-value argument. Note that the apply function automatically sets all image blocks outside the ROI to the default value of 0. To write the processed blocks to a TIFF file, specify the OutputLocation name-value argument with the target output file, and specify the Adapter name-value argument as a TIFF adapter.

outputFileMasked = "bostonparks_" + string(timestamp) + ".tif";
bimMasked = apply(bimCrop,@applyMask,ExtraImages=mask, ...
    blockLocationSet=bls,OutputLocation=outputFileMasked, ...
    Adapter=images.blocked.TIFF);

Display the mask in raster coordinates by using the bigimageshow function.

bigimageshow(bimMasked)

Convert TIFF File to GeoTIFF File

Convert the cropped and masked TIFF file to a valid GeoTIFF file by using the insertCroppedGeoTIFFTags helper function, which is attached to the example as a supporting file. The helper function generates and inserts the required GeoTIFF tags into the TIFF file. The data for the tags consists of coordinate information and additional metadata from the original GeoTIFF file.

insertCroppedGeoTIFFTags(outputFileMasked,"boston.tif",Rcrop)

To ensure that the cropped GeoTIFF file uses the correct map coordinates from the original GeoTIFF file, display the image from the cropped GeoTIFF file by using the mapshow function.

[bostonParks,RParks] = readgeoraster(outputFileMasked);
mapshow(bostonParks,RParks)

Helper Function

The applyMask helper function resizes a block of the mask to match the size of a block of image data, and sets the image data to 0 wherever the corresponding mask value is 0.

function imMasked = applyMask(data,mask)
    im = data.Data;
    mask = imresize(mask,size(im,[1 2]));
    mask = cast(mask,"like",im);
    imMasked = im.*mask;
end