Main Content

Find Water and Vegetation Regions in Spatially Referenced Multispectral Image

This example shows how to identify water and vegetation regions in a Landsat 8 multispectral image and spatially reference the image.

Spatial referencing is the process of assigning real-world coordinates to each pixel of a raster image to define its geospatial relationship. Mapping Toolbox™ enables you to obtain shapefile data within a geographic extent and display images on maps.

In this example, you use these steps to perform spatial referencing.

  1. Get the world coordinate limits.

  2. Get the size of the image.

  3. Create a reference raster cells to map coordinates using maprefcells (Mapping Toolbox) function.

Spectral indices characterize the specific features of interest of a target using its biophysical and chemical properties. These features of interest enable you to identify plant, water, and soil regions, as well as various forms of built-up regions. This example uses modified normalized difference water index (MNDWI) and green vegetation index (GVI) spectral indices to identify water and vegetation regions respectively. For more information on spectral indices, see Spectral Indices.

Load and Visualize Multispectral Data

Landsat 8 is an Earth observation satellite that carries the Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) instruments.

Download the Landsat 8 data set. The test data set has 8 spectral bands with wavelengths that range from 440 nm to 2200 nm. The test data is 7721-by-7651 pixels with a spatial resolution of 30 meters.

Download the data set and unzip the file by using the downloadLandsat8Dataset helper function. This function is attached to the example as a supporting file.

zipfile = "";
landsat8Data_url = "" + zipfile;
Downloading the Landsat 8 OLI dataset.
This can take several minutes to download...

Read the Landsat 8 multispectral data into the workspace as a hypercube object.

hCube = hypercube("LC08_L1TP_113082_20211206_20211206_01_RT_MTL.txt");

Estimate an RGB image from the data cube by using colorize function. Apply contrast stretching to enhance the contrast of the output RGB image.

rgbImg = colorize(hCube,Method="rgb",ContrastStretching=true);

Adjust image intensity values using the imadjustn function for better visualization.

rgbImg = imadjustn(rgbImg);

Display the RGB image of the test data. Notice that without spatial referencing, this figure does not provide any geographic information.

title("RGB Image of Data Cube")

Figure contains an axes object. The axes object with title RGB Image of Data Cube contains an object of type image.

Spatially Reference the Image

The Landsat 8 data set contains a GeoTIFF file. Obtain information about the GeoTIFF file by using the geotiffinfo (Mapping Toolbox) function.

filename = "LC08_L1TP_113082_20211206_20211206_01_RT_B1.TIF";
info = geotiffinfo(filename);

Calculate the latitude-longitude limits and world-xy limits using the corner coordinates of the GeoTIFF image.

cornerLat = info.CornerCoords.Lat;
cornerLon = info.CornerCoords.Lon;
cornerX = info.CornerCoords.X;
cornerY = info.CornerCoords.Y;

latLimits = [min(cornerLat(:)),max(cornerLat(:))];
lonLimits = [min(cornerLon(:)),max(cornerLon(:))];
xWorldLimits = [min(cornerX(:)),max(cornerX(:))];
yWorldLimits = [min(cornerY(:)),max(cornerY(:))];

Convert the GeoTIFF information to a map projection structure using geotiff2mstruct (Mapping Toolbox) function, to use for displaying the data in axesm-based map.

mstruct = geotiff2mstruct(info);

Import the shapefile worldcities.shp, which contains geographic information about major world cities as a geospatial table, by using the readgeotable (Mapping Toolbox) function. A geospatial table is a table or timetable object that contains a Shape variable and attribute variables. For more information on geospatial table, see Create Geospatial Tables (Mapping Toolbox). You can also use the worldrivers.shp and worldlakes.shp shapefiles to display major world rivers and lakes, respectively.

cities = readgeotable("worldcities.shp");

Create a polygon in geographic coordinates that represents the geographic extent of the region by using geopolyshape (Mapping Toolbox) object.

dataRegion = geopolyshape(latLimits([1 1 2 2 1]),lonLimits([1 2 2 1 1]));

Query the test data to determine which major cities are within the geographic extent of the rectangular data region. The data region contains a single city from the worldcities.shp geospatial table.

inRegion = isinterior(dataRegion,cities.Shape);
cityInd = find(inRegion);

Use basemaps, as well as with geographic axes and charts, to display the data region in axesm-based map, label the major city. You can specify different basemaps on which to plot the data using the geobasemap function. Use satellite as a basemap for this example.

basemap = "satellite";
hold on
for index = 1:cityInd
    text(cities(index,:).Shape.Latitude+0.07,cities(index,:).Shape.Longitude+0.03,cities(index,:).Name, ...
        HorizontalAlignment="left",FontWeight="bold",FontSize=14,Color=[1 1 1])
title("Geographic Extent of Landsat 8 Multispectal Image")

Figure contains an axes object. The axes object contains 231 objects of type polygon, scatter, text.

Find Water and Vegetation Regions in the Image

Compute the spectral index value for each pixel in the data cube by using spectralIndices function. Use the MNDWI and GVI to detect water and green vegetation regions, respectively.

indices = spectralIndices(hCube,["MNDWI","GVI"]);

Water regions typically have MNDWI values greater than 0. Vegetation regions typically have GVI values greater than 1. Specify the threshold values to perform thresholding of the MNDWI and GVI images to segment the water and green vegetation regions.

threshold = [0 1];

Generate binary images with a value of 1 for pixels with a score greater than the specified thresholds. All other pixels have a value of 0. The regions in the MNDWI and GVI binary images with a value of 1 correspond to the water and green vegetation regions, respectively.

Overlay the binary images on the RGB image by using labeloverlay function.

overlayImg = rgbImg;
labelColor = [0 0 1; 0 1 0];
for num = 1:numel(indices)
    indexMap = indices(num).IndexImage;
    thIndexMap = indexMap > threshold(num);
    overlayImg = labeloverlay(overlayImg,thIndexMap,Colormap=labelColor(num,:),Transparency=0.5);

Create a referencing object that specifies the world coordinate limits, and the raster size by using maprefcells (Mapping Toolbox) function.

R = maprefcells(xWorldLimits,yWorldLimits,size(overlayImg));

Resize the overlaid RGB image by using mapresize (Mapping Toolbox) function. For this example, reduce the size of overlaid RGB image to one-fourth of the original size.

scale = 1/4;
[reducedOverlayImg,reducedR] = mapresize(overlayImg,R,scale);

Display the overlaid image on an axesm-based map. The axes displays the water regions in blue and the green vegetation regions in green.

ax = axesm(mstruct,Grid="on", ...
    GColor=[1 1 1],GLineStyle="-", ...
    MapLatlimit=latLimits,MapLonLimit=lonLimits, ...
    ParallelLabel="on",PLabelLocation=0.5,PlabelMeridian="west", ...
    MeridianLabel="on",MlabelLocation=0.5,MLabelParallel="south", ...
    MLabelRound=-1,PLabelRound=-1, ...
    PLineVisible="on",PLineLocation=0.5, ...
axis off
dim = [0.8 0.5 0.3 0.3];
annotation(textbox=dim,String="Water Bodies", ...
    Color=[1 1 1], ...
    BackgroundColor=[0 0 1], ...
    FaceAlpha=0.5, ...
dim = [0.8 0.4 0.3 0.3];
annotation(textbox=dim,String="Green Vegetation", ...
    BackgroundColor=[0 1 0], ...
    FaceAlpha=0.5, ...
title("Water and Vegetation Region of Spatially Referenced Image")

Figure contains an axes object. The axes object with title Water and Vegetation Region of Spatially Referenced Image contains 12 objects of type image, line, text.

See Also

| | | (Mapping Toolbox) | (Mapping Toolbox) | (Mapping Toolbox)

Related Topics