Urban Form and Walkability, Kisumu, Kenya
Geospatial Data Science in Python

#| echo: true 

#| code-fold: true

# Stack all bands

def stack_bands(bands_list):

    stack = []

    for band_path in bands_list:

        with rio.open(band_path) as src:

            stack.append(src.read(1).flatten())  # Flatten each band

    return np.stack(stack, axis=1)

for band_path in bands_list:

    with rio.open(band_path) as src:

        array = src.read(1)

        print(array.shape)  # Print out the shape of each array

def stack_bands(bands_list):

    stack = []

    for band_path in bands_list:

        with rio.open(band_path) as src:

            array = src.read(1)

            # Only append bands that match the desired shape (279, 280) for example

            if array.shape == (279, 280):

                stack.append(array.flatten())  # Flatten each band

            else:

                print(f"Excluded {band_path} due to shape mismatch: {array.shape}")

    return np.vstack(stack).T  # Use vstack and transpose to get the correct shape for KMeans

# Get the stacked data excluding any mismatched bands

data_stack = stack_bands(bands_list)

# Stack all bands

def stack_bands(bands_list):

    stack = []

    for band_path in bands_list:

        with rio.open(band_path) as src:

            stack.append(src.read(1).flatten())  # Flatten each band

    return np.stack(stack, axis=1)

#| echo: true 

#| code-fold: true

# Reshape the data

data_stack = stack_bands(bands_list)

# Standardize the data

scaler = StandardScaler()

data_standardized = scaler.fit_transform(data_stack)

# Determine the optimal number of clusters

# Use methods like the Elbow method or Silhouette score here

# Sample for silhouette score

range_n_clusters = [2, 3, 4, 5, 6]

silhouette_avg = []

for num_clusters in range_n_clusters:

    kmeans = KMeans(n_clusters=num_clusters, random_state=42)

    cluster_labels = kmeans.fit_predict(data_standardized)

    silhouette_avg.append(silhouette_score(data_standardized, cluster_labels))

#| echo: true 

#| code-fold: true

# Choose the optimal number of clusters

optimal_clusters = range_n_clusters[silhouette_avg.index(max(silhouette_avg))]

# Apply K-Means Clustering

kmeans6 = KMeans(n_clusters=6, random_state=42)

y_pred = kmeans.fit_predict(data_standardized)

# Reshape classified data back to the original image shape

with rio.open(bands_list[0]) as src:

    meta = src.meta

    original_shape = (src.height, src.width)

classified_image6 = y_pred.reshape(original_shape)

# Visualize the result

plt.figure(figsize=(10, 10))

plt.imshow(classified_image6, cmap='terrain')

plt.colorbar()

plt.show()

# Assuming y_pred is your flattened array of predictions from KMeans

# and classified_image is your 2D array of cluster labels reshaped into the image's dimensions

# Count the occurrences of each cluster label

unique, counts = np.unique(classified_image6, return_counts=True)

cluster_counts6 = dict(zip(unique, counts))

# Calculate the total number of pixels

total_pixels6 = classified_image6.size

# Calculate the percentage of each cluster

cluster_percentages6 = {k: (v / total_pixels6) * 100 for k, v in cluster_counts6.items()}

# Print the percentages

for cluster, percentage in cluster_percentages6.items():

    print(f'Cluster {cluster}: {percentage:.2f}%')

# Apply K-Means Clustering

kmeans3 = KMeans(n_clusters=3, random_state=42)

y_pred = kmeans3.fit_predict(data_standardized)

# Reshape classified data back to the original image shape

with rio.open(bands_list[0]) as src:

    meta = src.meta

    original_shape = (src.height, src.width)

classified_image3 = y_pred.reshape(original_shape)

# Visualize the result

plt.figure(figsize=(10, 10))

plt.imshow(classified_image3, cmap='terrain')

plt.colorbar()

plt.show()

# Assuming y_pred is your flattened array of predictions from KMeans

# and classified_image is your 2D array of cluster labels reshaped into the image's dimensions

# Count the occurrences of each cluster label

unique, counts = np.unique(classified_image3, return_counts=True)

cluster_counts3 = dict(zip(unique, counts))

# Calculate the total number of pixels

total_pixels3 = classified_image3.size

# Calculate the percentage of each cluster

cluster_percentages3 = {k: (v / total_pixels3) * 100 for k, v in cluster_counts3.items()}

# Print the percentages

for cluster, percentage in cluster_percentages3.items():

    print(f'Cluster {cluster}: {percentage:.2f}%')

clustered_models = ClusteredBands(complete_dataset)

clustered_models.set_raster_stack()