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.
Get the world coordinate limits.
Get the size of the image.
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 = "LC08_L1TP_113082_20211206_20211206_01_RT.zip"; landsat8Data_url = "https://ssd.mathworks.com/supportfiles/image/data/" + zipfile; downloadLandsat8Dataset(landsat8Data_url,pwd);
Downloading the Landsat 8 OLI dataset. This can take several minutes to download... Done.
Read the Landsat 8 multispectral data into the workspace as a
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.
figure imshow(rgbImg) title("RGB Image of Data Cube")
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
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
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
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.
figure basemap = "satellite"; geobasemap(basemap) hold on geoplot(dataRegion,LineWidth=2,EdgeColor="yellow",FaceColor="red",FaceAlpha=0.1); geoscatter(cities(cityInd,:).Shape.Latitude,cities(cityInd,:).Shape.Longitude,"filled") 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]) end title("Geographic Extent of Landsat 8 Multispectal Image")
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
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); end
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.
figure 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, ... MLineVisible="on",MlineLocation=0.5); mapshow(reducedOverlayImg,reducedR) 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, ... FitBoxToText="on"); dim = [0.8 0.4 0.3 0.3]; annotation(textbox=dim,String="Green Vegetation", ... BackgroundColor=[0 1 0], ... FaceAlpha=0.5, ... FitBoxToText="on"); title("Water and Vegetation Region of Spatially Referenced Image")