05. Calculating landcover change over time

imports

import leafmap
import leafmap.colormaps as cm
import requests
import rasterio as rio
from rasterio.merge import merge
import glob
from rasterio.plot import reshape_as_image
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

searching for data

url = "https://earth-search.aws.element84.com/v1/"
collection = "sentinel-2-l2a"
time_range = "2023-08-01/2023-08-31"

# bbox for dallas metro
bbox = [
    -97.06213756027009,
    32.97324551867027,
    -96.46807822577594,
    33.3578329610085,
]

search_gdf = leafmap.stac_search(
    url=url,
    max_items=10,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    query={"eo:cloud_cover": {"lt": 20}},
    sortby=[{"field": "properties.eo:cloud_cover", "direction": "asc"}],
    get_gdf=True,
)

search_gdf.columns

Index(['geometry', 'created', 'platform', 'constellation', 'instruments',
       'eo:cloud_cover', 'mgrs:utm_zone', 'mgrs:latitude_band',
       'mgrs:grid_square', 'grid:code', 'view:sun_azimuth',
       'view:sun_elevation', 's2:degraded_msi_data_percentage',
       's2:nodata_pixel_percentage', 's2:saturated_defective_pixel_percentage',
       's2:dark_features_percentage', 's2:cloud_shadow_percentage',
       's2:vegetation_percentage', 's2:not_vegetated_percentage',
       's2:water_percentage', 's2:unclassified_percentage',
       's2:medium_proba_clouds_percentage', 's2:high_proba_clouds_percentage',
       's2:thin_cirrus_percentage', 's2:snow_ice_percentage',
       's2:product_type', 's2:processing_baseline', 's2:product_uri',
       's2:generation_time', 's2:datatake_id', 's2:datatake_type',
       's2:datastrip_id', 's2:granule_id', 's2:reflectance_conversion_factor',
       'datetime', 's2:sequence', 'earthsearch:s3_path',
       'earthsearch:payload_id', 'earthsearch:boa_offset_applied',
       'processing:software', 'updated', 'proj:code'],
      dtype='str')

search_gdf

exploring and filtering search results

search_gdf.plot("mgrs:grid_square", alpha=0.25)

<Axes: >

search_gdf["s2:granule_id"].value_counts()

s2:granule_id
S2B_OPER_MSI_L2A_TL_2BPS_20230819T223356_A033700_T14SQB_N05.09    1
S2A_OPER_MSI_L2A_TL_2APS_20230804T235402_A042394_T14SQB_N05.09    1
S2A_OPER_MSI_L2A_TL_2APS_20230818T002256_A042580_T14SPC_N05.09    1
S2B_OPER_MSI_L2A_TL_2BPS_20230819T223356_A033700_T14SPB_N05.09    1
S2A_OPER_MSI_L2A_TL_2APS_20230815T002553_A042537_T14SPB_N05.09    1
S2A_OPER_MSI_L2A_TL_2APS_20230815T002553_A042537_T14SQB_N05.09    1
S2B_OPER_MSI_L2A_TL_2BPS_20230819T223356_A033700_T14SPC_N05.09    1
S2A_OPER_MSI_L2A_TL_2APS_20230818T002256_A042580_T14SQC_N05.09    1
S2A_OPER_MSI_L2A_TL_2APS_20230804T235402_A042394_T14SPB_N05.09    1
S2B_OPER_MSI_L2A_TL_2BPS_20230819T223356_A033700_T14SQC_N05.09    1
Name: count, dtype: int64

search_gdf["area"] = search_gdf.geometry.area

search_gdf.sort_values("area", ascending=False, inplace=True)

search_gdf[["area"]]

search_gdf.drop_duplicates("mgrs:grid_square", inplace=True)

search_gdf.iloc[0:4].plot("s2:granule_id", alpha=0.5)

<Axes: >

# # create a query searching for unique datastrip ids based on the above
q = {"s2:granule_id": {"in": search_gdf["s2:granule_id"].iloc[0:4].unique().tolist()}}

search_gdf2 = leafmap.stac_search(
    url=url,
    max_items=4,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    query=q,
    get_gdf=True,
)

search_gdf2

search_gdf2.plot("mgrs:grid_square", alpha=0.5)

<Axes: >

plot to accentuate paved and non-paved divide

search_gdf2_list = leafmap.stac_search(
    url=url,
    max_items=4,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    query=q,
    get_links=True,
)

search_gdf2_list

['https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2B_14SQB_20230819_0_L2A',
 'https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2B_14SQC_20230819_0_L2A',
 'https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2A_14SPC_20230817_0_L2A',
 'https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2A_14SPB_20230814_0_L2A']

download files for analysis

def get_item_json(item_url: str) -> dict:
    r = requests.get(item_url, timeout=30)
    r.raise_for_status()
    return r.json()


def get_raster_band_urls(item_url: str, bands: list | None = None) -> dict:
    """Return {band_name: cog_url} for a STAC item, no titiler dependency."""
    stac = get_item_json(item_url)
    band_urls = {
        key: asset["href"]
        for key, asset in stac["assets"].items()
        if asset["href"].startswith("http")
        and "image/tiff"
        in asset.get("type", "")  # COGs only; skips thumbnail jpg, xml metadata, jp2
    }
    if bands:
        band_urls = {b: band_urls[b] for b in bands if b in band_urls}
    return band_urls


def download_stac_layers(layers, out_dir, bands=None):
    for layer in layers:
        band_urls = get_raster_band_urls(layer, bands)
        if not band_urls:
            print(f"WARNING: no COG assets found for {layer}")
            continue
        for band, url in band_urls.items():
            out_file = f"{out_dir}/{layer.split('/')[-1]}_{band}.tif"
            leafmap.download_file(url, out_file, overwrite=False)


def get_stac_crs(item_url: str) -> int:
    """Handle both proj:epsg (projection ext v1) and proj:code (v2, e.g. 'EPSG:32614')."""
    stac = get_item_json(item_url)
    props = stac.get("properties", {})
    if "proj:epsg" in props:
        return props["proj:epsg"]
    if "proj:code" in props:
        return int(props["proj:code"].split(":")[-1])
    # some catalogs moved projection info to the asset level
    for asset in stac["assets"].values():
        if "proj:code" in asset:
            return int(asset["proj:code"].split(":")[-1])
        if "proj:epsg" in asset:
            return asset["proj:epsg"]
    raise KeyError(f"No projection metadata found on {item_url}")

stac_crs = get_stac_crs(search_gdf2_list[0])

stac_crs

32614

download_stac_layers(
    search_gdf2_list,
    "../Data/stac/dallas",
)

../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_aot.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_blue.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_coastal.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_green.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_nir.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_nir08.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_nir09.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_red.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_rededge1.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_rededge2.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_rededge3.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_scl.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_swir16.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_swir22.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_visual.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQB_20230819_0_L2A_wvp.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_aot.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_blue.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_coastal.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_green.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_nir.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_nir08.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_nir09.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_red.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_rededge1.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_rededge2.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_rededge3.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_scl.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_swir16.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_swir22.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_visual.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2B_14SQC_20230819_0_L2A_wvp.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_aot.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_blue.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_coastal.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_green.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_nir.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_nir08.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_nir09.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_red.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_rededge1.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_rededge2.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_rededge3.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_scl.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_swir16.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_swir22.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_visual.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPC_20230817_0_L2A_wvp.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_aot.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_blue.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_coastal.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_green.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_nir.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_nir08.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_nir09.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_red.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_rededge1.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_rededge2.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_rededge3.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_scl.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_swir16.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_swir22.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_visual.tif already exists. Skip downloading. Set overwrite=True to overwrite.
../Data/stac/dallas/S2A_14SPB_20230814_0_L2A_wvp.tif already exists. Skip downloading. Set overwrite=True to overwrite.

# mosaic files based on band name


def mosaic_by_band(
    dir,
    bands: list = [
        "red",
        "blue",
        "green",
        "nir",
        "coastal",
        "nir08",
        "nir09",
        "rededge1",
        "rededge2",
        "rededge3",
        "scl",
        "swir16",
        "swir22",
    ],
    crs: str = None,
):
    mosaics = {}
    for band in bands:
        files = glob.glob(f"{dir}/*{band}.tif")

        out_file = f"{dir}/mosaic_{band}_.tif"
        raster_data = [rio.open(f) for f in files]

        mosaic, out_trans = merge(raster_data)
        out_meta = raster_data[0].meta.copy()
        out_meta.update(
            {
                "driver": "GTiff",
                "height": mosaic.shape[1],
                "width": mosaic.shape[2],
                "transform": out_trans,
                "crs": f"epsg:{crs}",
            }
        )

        with rio.open(out_file, "w", **out_meta) as dest:
            dest.write(mosaic)

        mosaics[band] = out_file

    return mosaics

mosaic_bands_all = mosaic_by_band("../Data/stac/dallas", crs=stac_crs)

mosaic_bands_all

{'red': '../Data/stac/dallas/mosaic_red_.tif',
 'blue': '../Data/stac/dallas/mosaic_blue_.tif',
 'green': '../Data/stac/dallas/mosaic_green_.tif',
 'nir': '../Data/stac/dallas/mosaic_nir_.tif',
 'coastal': '../Data/stac/dallas/mosaic_coastal_.tif',
 'nir08': '../Data/stac/dallas/mosaic_nir08_.tif',
 'nir09': '../Data/stac/dallas/mosaic_nir09_.tif',
 'rededge1': '../Data/stac/dallas/mosaic_rededge1_.tif',
 'rededge2': '../Data/stac/dallas/mosaic_rededge2_.tif',
 'rededge3': '../Data/stac/dallas/mosaic_rededge3_.tif',
 'scl': '../Data/stac/dallas/mosaic_scl_.tif',
 'swir16': '../Data/stac/dallas/mosaic_swir16_.tif',
 'swir22': '../Data/stac/dallas/mosaic_swir22_.tif'}

nir = rio.open(mosaic_bands_all["nir"]).read(1).astype("float32")
red = rio.open(mosaic_bands_all["red"]).read(1).astype("float32")

calculate NDVI

ndvi = (nir - red) / (nir + red)

# fill nan with -2, outside of the range but not so far outside
ndvi[np.isnan(ndvi)] = -1.1

ndvi

array([[ 0.36826923,  0.36826923,  0.36826923, ...,  0.72460496,
         0.7580539 ,  0.7661823 ],
       [ 0.36826923,  0.36826923,  0.36826923, ...,  0.75057286,
         0.75632185,  0.7632707 ],
       [ 0.36826923,  0.36826923,  0.36826923, ...,  0.75344265,
         0.7732096 ,  0.7714475 ],
       ...,
       [-1.1       , -1.1       , -1.1       , ...,  0.28488854,
         0.29016393,  0.29233912],
       [-1.1       , -1.1       , -1.1       , ...,  0.2907453 ,
         0.2860676 ,  0.29248366],
       [-1.1       , -1.1       , -1.1       , ...,  0.29306123,
         0.27983367,  0.277882  ]], shape=(20982, 20976), dtype=float32)

ndvi_mult = ndvi * 10  # to be able to visualize more easily

ndvi_mult

array([[  3.6826923,   3.6826923,   3.6826923, ...,   7.24605  ,
          7.5805387,   7.6618233],
       [  3.6826923,   3.6826923,   3.6826923, ...,   7.5057287,
          7.5632186,   7.6327066],
       [  3.6826923,   3.6826923,   3.6826923, ...,   7.5344267,
          7.7320957,   7.7144747],
       ...,
       [-11.       , -11.       , -11.       , ...,   2.8488853,
          2.9016395,   2.923391 ],
       [-11.       , -11.       , -11.       , ...,   2.9074528,
          2.860676 ,   2.9248366],
       [-11.       , -11.       , -11.       , ...,   2.9306123,
          2.7983367,   2.77882  ]], shape=(20982, 20976), dtype=float32)

ndvi_image = leafmap.array_to_image(ndvi, source=mosaic_bands_all["nir"])

ndvi_multi_image = leafmap.array_to_image(ndvi_mult, source=mosaic_bands_all["nir"])

m = leafmap.Map()
m.add_raster(ndvi_multi_image, layer_name="NDVI", vmin=-10, vmax=10, cmap="PiYG")
m

Map(center=[33.371449, -96.78545], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title',…

# get histogram of ndvi image values
plt.hist(ndvi[ndvi != -1.1].ravel(), bins=100, color="red", alpha=0.7)

(array([1.0000000e+00, 1.0000000e+00, 0.0000000e+00, 0.0000000e+00,
        0.0000000e+00, 1.0000000e+00, 0.0000000e+00, 2.0000000e+00,
        5.0000000e+00, 1.8000000e+01, 2.7000000e+01, 9.5000000e+01,
        2.7300000e+02, 5.7100000e+02, 9.5100000e+02, 1.4730000e+03,
        2.2680000e+03, 3.1480000e+03, 4.4330000e+03, 6.0950000e+03,
        8.7460000e+03, 1.2123000e+04, 1.5750000e+04, 1.9139000e+04,
        2.3716000e+04, 2.8782000e+04, 3.4653000e+04, 4.5633000e+04,
        6.4668000e+04, 7.6474000e+04, 8.9700000e+04, 1.2927000e+05,
        1.6543500e+05, 2.3599500e+05, 3.1663400e+05, 3.8703400e+05,
        4.7900100e+05, 6.2561400e+05, 7.5519000e+05, 9.0728100e+05,
        1.0850600e+06, 1.1202380e+06, 8.8427400e+05, 6.4577800e+05,
        6.0382200e+05, 5.3926600e+05, 3.4452000e+05, 3.3818100e+05,
        3.4321800e+05, 4.3010100e+05, 5.9539900e+05, 7.7239800e+05,
        1.0676070e+06, 1.3586870e+06, 1.6881750e+06, 2.1018380e+06,
        2.6668380e+06, 3.6052190e+06, 4.7291470e+06, 6.0544950e+06,
        7.6989430e+06, 9.4903910e+06, 1.1130270e+07, 1.2534731e+07,
        1.3674737e+07, 1.4563089e+07, 1.5116707e+07, 1.5522868e+07,
        1.5675501e+07, 1.5614456e+07, 1.5352237e+07, 1.4917200e+07,
        1.4369598e+07, 1.3655650e+07, 1.2803444e+07, 1.1886788e+07,
        1.0984785e+07, 1.0146013e+07, 9.3278850e+06, 8.6188730e+06,
        7.9977950e+06, 7.4678020e+06, 7.0248890e+06, 6.6651540e+06,
        6.4014320e+06, 6.2831160e+06, 6.3451950e+06, 6.6194690e+06,
        7.1241650e+06, 7.7236790e+06, 7.7266050e+06, 5.2941400e+06,
        1.6410080e+06, 1.6154700e+05, 8.3670000e+03, 7.4200000e+02,
        1.5000000e+01, 8.0000000e+00, 5.0000000e+00, 2.6400000e+02]),
 array([-0.99696511, -0.97700125, -0.95703739, -0.93707347, -0.91710961,
        -0.89714575, -0.87718189, -0.85721803, -0.83725417, -0.81729025,
        -0.79732639, -0.77736253, -0.75739866, -0.73743474, -0.71747088,
        -0.69750702, -0.67754316, -0.6575793 , -0.63761544, -0.61765158,
        -0.59768766, -0.5777238 , -0.55775994, -0.53779602, -0.51783216,
        -0.49786833, -0.47790444, -0.45794058, -0.43797672, -0.41801286,
        -0.39804894, -0.37808508, -0.35812122, -0.33815736, -0.3181935 ,
        -0.29822958, -0.27826571, -0.25830185, -0.23833799, -0.21837413,
        -0.19841021, -0.17844635, -0.15848249, -0.13851863, -0.11855477,
        -0.09859091, -0.07862699, -0.05866313, -0.03869927, -0.01873541,
         0.00122845,  0.02119237,  0.04115623,  0.06112009,  0.08108395,
         0.10104781,  0.12101167,  0.14097553,  0.1609394 ,  0.18090338,
         0.20086724,  0.2208311 ,  0.24079496,  0.26075882,  0.28072268,
         0.30068654,  0.3206504 ,  0.34061426,  0.36057812,  0.38054198,
         0.40050596,  0.42046982,  0.44043368,  0.46039754,  0.4803614 ,
         0.50032526,  0.52028912,  0.54025298,  0.56021684,  0.5801807 ,
         0.60014468,  0.62010854,  0.64007241,  0.66003627,  0.68000013,
         0.69996399,  0.71992785,  0.73989171,  0.75985557,  0.77981943,
         0.79978329,  0.81974727,  0.83971113,  0.85967499,  0.87963885,
         0.89960271,  0.91956657,  0.93953043,  0.95949429,  0.97945815,
         0.99942213]),
 <BarContainer object of 100 artists>)

def categorize_np_array(array):

    cat_array = np.zeros(array.shape)
    cat_array[array < 1] = 5
    cat_array[array < 0.8] = 4
    cat_array[array < 0.6] = 3
    cat_array[array < 0.4] = 2
    cat_array[array < 0.2] = 1

    return cat_array

cat_ndvi = categorize_np_array(ndvi)

cat_ndvi

array([[2., 2., 2., ..., 4., 4., 4.],
       [2., 2., 2., ..., 4., 4., 4.],
       [2., 2., 2., ..., 4., 4., 4.],
       ...,
       [1., 1., 1., ..., 2., 2., 2.],
       [1., 1., 1., ..., 2., 2., 2.],
       [1., 1., 1., ..., 2., 2., 2.]], shape=(20982, 20976))

ndvi_image_cat = leafmap.array_to_image(cat_ndvi, source=mosaic_bands_all["nir"])

m = leafmap.Map()
m.add_raster(
    ndvi_image_cat, layer_name="NDVI", vmin=1, vmax=5, cmap="Set1", draw_control=False
)
m.add_colormap(
    "Set1",
    label="ndvi",
    width=8.0,
    height=0.4,
    orientation="horizontal",
    vmin=1,
    vmax=5,
)

m

Map(center=[33.371449, -96.78545], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title',…

k-means

def resample_raster(in_file, out_file, upscale_factor=0.5):
    with rio.open(in_file) as dataset:

        # resample data to target shape
        data = dataset.read(
            out_shape=(
                dataset.count,
                int(dataset.height * upscale_factor),
                int(dataset.width * upscale_factor),
            ),
            resampling=rio.enums.Resampling.bilinear,
        )

        # scale image transform
        transform = dataset.transform * dataset.transform.scale(
            (dataset.width / data.shape[-1]), (dataset.height / data.shape[-2])
        )

        # update metadata
        meta = dataset.meta.copy()
        meta.update(
            {
                "height": data.shape[1],
                "width": data.shape[2],
                "transform": transform,
            }
        )

        with rio.open(out_file, "w", **meta) as dst:
            dst.write(data)

for band in mosaic_bands_all.keys():
    print(band, rio.open(mosaic_bands_all[band]).read(1).shape)

red (20982, 20976)
blue (20982, 20976)
green (20982, 20976)
nir (20982, 20976)
coastal (3497, 3496)
nir08 (10491, 10488)
nir09 (3497, 3496)
rededge1 (10491, 10488)
rededge2 (10491, 10488)
rededge3 (10491, 10488)
scl (10491, 10488)
swir16 (10491, 10488)
swir22 (10491, 10488)

# set resampling ratio for each band- note that some are higher resolution than others and need a different scaling factor
bands_to_resample = {
    "red": 0.25,
    "blue": 0.25,
    "green": 0.25,
    "nir": 0.25,
    "nir08": 0.5,
    "rededge1": 0.5,
    "rededge2": 0.5,
    "rededge3": 0.5,
    "scl": 0.5,
    "swir16": 0.5,
    "swir22": 0.5,
}

# resample the bands so they are all smaller and uniform
for band, ratio in bands_to_resample.items():
    resample_raster(
        mosaic_bands_all[band], f"../Data/stac/dallas/upscaled_{band}.tif", ratio
    )
    mosaic_bands_all[f"{band}"] = f"../Data/stac/dallas/upscaled_{band}.tif"

mosaic_bands_all

{'red': '../Data/stac/dallas/upscaled_red.tif',
 'blue': '../Data/stac/dallas/upscaled_blue.tif',
 'green': '../Data/stac/dallas/upscaled_green.tif',
 'nir': '../Data/stac/dallas/upscaled_nir.tif',
 'coastal': '../Data/stac/dallas/mosaic_coastal_.tif',
 'nir08': '../Data/stac/dallas/upscaled_nir08.tif',
 'nir09': '../Data/stac/dallas/mosaic_nir09_.tif',
 'rededge1': '../Data/stac/dallas/upscaled_rededge1.tif',
 'rededge2': '../Data/stac/dallas/upscaled_rededge2.tif',
 'rededge3': '../Data/stac/dallas/upscaled_rededge3.tif',
 'scl': '../Data/stac/dallas/upscaled_scl.tif',
 'swir16': '../Data/stac/dallas/upscaled_swir16.tif',
 'swir22': '../Data/stac/dallas/upscaled_swir22.tif'}

# create a new dict with all bands excpct coastal
mosaic_bands_all_no_coastal = {
    k: v for k, v in mosaic_bands_all.items() if k not in ["coastal", "nir09"]
}

mosaic_bands_all_no_coastal

{'red': '../Data/stac/dallas/upscaled_red.tif',
 'blue': '../Data/stac/dallas/upscaled_blue.tif',
 'green': '../Data/stac/dallas/upscaled_green.tif',
 'nir': '../Data/stac/dallas/upscaled_nir.tif',
 'nir08': '../Data/stac/dallas/upscaled_nir08.tif',
 'rededge1': '../Data/stac/dallas/upscaled_rededge1.tif',
 'rededge2': '../Data/stac/dallas/upscaled_rededge2.tif',
 'rededge3': '../Data/stac/dallas/upscaled_rededge3.tif',
 'scl': '../Data/stac/dallas/upscaled_scl.tif',
 'swir16': '../Data/stac/dallas/upscaled_swir16.tif',
 'swir22': '../Data/stac/dallas/upscaled_swir22.tif'}

mosaic_bands_all_no_coastal_cols = mosaic_bands_all_no_coastal.keys()

for band in mosaic_bands_all_no_coastal.keys():
    print(band, rio.open(mosaic_bands_all_no_coastal[band]).read(1).shape)

mosaic_bands_all_no_coastal = np.stack(
    [
        rio.open(mosaic_bands_all_no_coastal[band]).read(1)
        for band in mosaic_bands_all_no_coastal.keys()
    ],
    axis=-1,
)

red (5245, 5244)
blue (5245, 5244)
green (5245, 5244)
nir (5245, 5244)
nir08 (5245, 5244)
rededge1 (5245, 5244)
rededge2 (5245, 5244)
rededge3 (5245, 5244)
scl (5245, 5244)
swir16 (5245, 5244)
swir22 (5245, 5244)

reshaped_img = reshape_as_image(mosaic_bands_all_no_coastal)

reshaped_img.shape

(5244, 11, 5245)

mosaic_bands_all_no_coastal.shape

(5245, 5244, 11)

rows, cols, bands = mosaic_bands_all_no_coastal.shape

k = 10  # num of clusters

kmeans_predictions = KMeans(n_clusters=k, random_state=0).fit(
    mosaic_bands_all_no_coastal.reshape(rows * cols, bands)
)

kmeans_predictions_2d = kmeans_predictions.labels_.reshape(rows, cols)

kmeans_predictions_2d

array([[7, 3, 4, ..., 5, 5, 5],
       [7, 4, 3, ..., 5, 5, 5],
       [9, 3, 3, ..., 5, 5, 5],
       ...,
       [1, 1, 1, ..., 4, 4, 4],
       [1, 1, 1, ..., 4, 4, 3],
       [1, 1, 1, ..., 4, 3, 3]], shape=(5245, 5244), dtype=int32)

# plot prediction results

m = leafmap.Map()
m.add_raster(
    leafmap.array_to_image(
        kmeans_predictions_2d,
        source=mosaic_bands_all["red"],
    ),
    layer_name="KMeans",
    colormap="Set1",
    draw_control=False,
    vmin=0,
    vmax=9,
)

# add a legend
m.add_colormap(
    "Set1",
    label="KMeans",
    width=8.0,
    height=0.4,
    orientation="horizontal",
    vmin=0,
    vmax=9,
)

m

Map(center=[33.371449, -96.78545], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title',…

xb, yb = "red", "nir"  # two test bands
labels = kmeans_predictions_2d.reshape(-1)

df = pd.DataFrame(
    mosaic_bands_all_no_coastal.reshape(rows * cols, bands)[idx],
    columns=mosaic_bands_all_no_coastal_cols,
)
df["cluster"] = labels[idx]
print(df.columns.to_list())
g = sns.FacetGrid(
    df,
    col="cluster",
    col_wrap=4,
    height=2.6,
    aspect=1,
)


def panel(data, color, **kws):
    ax = plt.gca()
    # reference points
    ax.scatter(
        df[xb],
        df[yb],
        s=2,
        color="0.85",
        linewidth=0,
        rasterized=True,
    )

    k = data["cluster"].iloc[0]
    ax.scatter(data[xb], data[yb], s=3, color=palette[k], alpha=0.5, linewidth=0)


g.map_dataframe(panel)
g.set_axis_labels(f"{xb} (DN)", f"{yb} (DN)")
g.set_titles("cluster {col_name}")
g.figure.suptitle("NIR vs Red by cluster", y=1.03)

['red', 'blue', 'green', 'nir', 'nir08', 'rededge1', 'rededge2', 'rededge3', 'scl', 'swir16', 'swir22', 'cluster']

Text(0.5, 1.03, 'NIR vs Red by cluster')