Assign7 Classification problem – solution

67. Assign7 Classification problem – solution#

Make a jupyter notebook that reproduces the false color examples inLandsat 8 false color examples for your scene (you’ll need to rerun Clipping multiple bands– v0.3 version 3 to get the clipped fmask, and rerun Landsat 1: Dowloading Landsat and Sentinel data from NASA to download all of the HLS tifs for bands 1,2,3,4,5,6,7,9,10,11,fmask if you don’t have them).
Choose one band combination that looks interesting, and compare it with the land classification you created using Working with surface class data for your image with the same bounding box and pixel size – comment on any similarities and differences you can find. Is the classification accurate?

import xarray
import rioxarray
from matplotlib import pyplot as plt
import numpy as np
from skimage import exposure, img_as_ubyte
from IPython.display import Image
from pathlib import Path
from numpy.typing import NDArray
from rasterio.enums import Resampling
import seaborn as sns

67.1. False color bands for Vancouver#

bands={'B01':'Coastal_Aerosol',
        'B02':'Blue',
        'B03':'Green',
        'B04':'Red',
        'B05':'NIR',
        'B06':'SWIR1',
        'B07':'SWIR2',
        'B09':'Cirrus',
        'B10':'TIRS1',
        'B11':'TIRS2',
        'fmask':'fmask'}

67.1.1. Read all bands into a dictionary#

Use the bands dictionary to identify the band by its name (Blue, Green, etc) and store it in scene_dict. Masking fmask would convert it from 8 bit to float, so we need to special-case the fmask file.

data_dir = Path().home() / 'repos/a301/satdata/landsat'
all_tifs = list(data_dir.glob('**/week10*clipped_*.tif'))
scene_dict = {}
for key,bandname in bands.items():
    band_tif = None
    for the_tif in all_tifs:
        if str(the_tif).find(bandname) > -1:
            print(f"reading {key}:{the_tif}")
            if key == 'fmask':
                scene_dict[key] = rioxarray.open_rasterio(the_tif)
            else:
                scene_dict[key] = rioxarray.open_rasterio(the_tif, mask_and_scale=True)
            continue
   
    

reading B01:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_Coastal_Aerosol.tif
reading B02:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_Blue.tif
reading B03:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_Green.tif
reading B04:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_Red.tif
reading B05:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_NIR.tif
reading B06:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_SWIR1.tif
reading B07:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_SWIR2.tif
reading B09:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_Cirrus.tif
reading B10:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_TIRS1.tif
reading B11:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_TIRS2.tif
reading fmask:/Users/phil/repos/a301/satdata/landsat/vancouver_2023/week10/week10_clipped_fmask.tif

67.1.3. Make the mask#

bool_image = make_bool_mask(scene_dict['fmask'])
plt.imshow(bool_image.squeeze())

<matplotlib.image.AxesImage at 0x13f962900>

../_images/2bbfd646cad5adeef4de476e9fa49e6f2203a25bf1faf280be1c8825cec92850.png

67.2. Urban swir-2, swir-1, red#

concrete and bare soil have approximately constant reflectivities between 1.6 and 2.2 microns, while vegetation reflects more in swir-1 than swir-2. This combination distinguishes between types of urban development. Less contrast for turbid/fresh water. See: band764 detail

urban = make_false_color(ds_allbands, band_names=["B07","B06","B04"])
fig5, ax5 = plt.subplots(1,1,figsize=(6,9))
urban.plot.imshow(ax=ax5);
ax5.set(title="Urban 764");

../_images/accacacd09935f70cd5d4a486bd5ede7d65115a70c7fb4cf7118d5db371887e5.png

67.3. Read in the classification and resample#

The land use/land classification process assigns a category between 0-11 for the following 9 categories (they have combined former categories 3 and 6 in the the category 11 (rangeland).

\[\begin{split} \begin{array}{|ll|} \hline 0 & \text { No data } \\ \hline 1 & \text { Water } \\ \hline 2 & \text { Trees } \\ \hline 4 & \text { Flooded Vegetation } \\ \hline 5 & \text { Crops } \\ \hline 7 & \text { Built Area } \\ \hline 8 & \text { Bare Ground } \\ \hline 9 & \text { Snow/lce } \\ \hline 10 & \text { Clouds } \\ \hline 11 & \text { Rangeland } \\ \hline \end{array} \end{split}\]

67.3.1. get the bounds from band7#

the_tif = list(data_dir.glob('**/vancouver_2023/10U*2023*tif'))[0]
print(the_tif)

/Users/phil/repos/a301/satdata/landsat/vancouver_2023/10U_2023.tif

67.3.2. get the classification file#

band7 = ds_allbands['B07'].squeeze()
bounds = band7.rio.bounds()
land_class = rioxarray.open_rasterio(the_tif)
land_class = land_class.squeeze()

67.4. resample to band 7 grid#

matched_class = land_class.rio.reproject_match(band7,
                                                  resampling=Resampling.mode)

matched_class.plot.imshow();

../_images/5c12f6210aa722f39cbb866c6822483a03168a20d595ce002234edade7278e8d.png

matched_class.plot.hist()

(array([1.82188e+05, 1.23770e+04, 0.00000e+00, 1.66000e+02, 9.64000e+02,
        0.00000e+00, 3.63020e+04, 1.31000e+02, 0.00000e+00, 7.87200e+03]),
 array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11.]),
 <BarContainer object of 10 artists>)

../_images/09985dcb0c39a9301149dacabd23f2f93f826ce6e0978f9add88678cfc3c1a2e.png

67.5. Discussion#

The classifier appeaars to be doing a good job of distinguishing the concrete runways as built, and the mowed grass around the runways as rangeland. for the YVR runways. It also picks out the UBC golf course as rangeland, but seems to also classify some sandy areas as rangeland around Iona spit. Not sure what’s going on with the crop classification in Point Grey. Would need to look at sentinel msi bands which are used to develop the classifier to see if they are making some other kind of distinction in their bands 5-7 which landsat doesn’t have

hit_built = matched_class.data == 7
hit_range = matched_class.data ==11
band6 = ds_allbands['B06'].squeeze()
band4 = ds_allbands['B04'].squeeze()

band7.shape

(600, 400)

band7_built = band7.data[hit_built]
band6_built = band6.data[hit_built]
band4_built = band4.data[hit_built]
band7_range = band7.data[hit_range]
band6_range = band6.data[hit_range]
band4_range = band4.data[hit_range]

67.5.1. Jointplot for built pixels bands 6,7#

fig = sns.jointplot(
    x=band6_built,
    y=band7_built,
    kind="hex",
    xlim=(0.08, 0.3),
    ylim=(0.0, 0.24),
    color="#4CB391",
    gridsize=100,
)
fig.set_axis_labels("band 6 built", "band 7 built");

../_images/34dd40f61eeb25103cd2c74ddd1fbf68d595a712376f8411f0b0d0088f4c15b2.png

67.5.2. Joint plot for range pixels, bands 6 and 7#

fig = sns.jointplot(
    x=band6_range,
    y=band7_range,
    kind="hex",
    xlim=(0.08, 0.3),
    ylim=(0.0, 0.24),
    color="#4CB391",
    gridsize=100,
)
fig.set_axis_labels("band 6 range", "band 7 range");

../_images/edc5b0025f63059b3e60a2d1c041e2a8e84c0efb7d6819310f45aad18599fdd7.png

67.5.3. Jointplot for built pixels bands 4,7#

fig = sns.jointplot(
    x=band4_built,
    y=band7_built,
    kind="hex",
    xlim=(0.0, 0.15),
    ylim=(0.0, 0.24),
    color="#4CB391",
    gridsize=100,
)
fig.set_axis_labels("band 4 range", "band 7 range");

../_images/1a147a450ef3d8434345603cb5b2a6b7e57fc158e688874d388ae93ce3f411c8.png

67.5.4. Joint plot for range pixels bands 4, 7#

fig = sns.jointplot(
    x=band4_range,
    y=band7_range,
    kind="hex",
    xlim=(0.0, 0.15),
    ylim=(0.0, 0.24),
    color="#4CB391",
    gridsize=100,
)
fig.set_axis_labels("band 4 range", "band 7 range");

../_images/df259d22b512ba7bdfd625110d343c2ba8bb370402c5858a7e60a9002061c950.png