Forest & Buildings, Removed DEM Error

import ee
import colorcet
from matplotlib.cm import get_cmap
import geemap
from datetime import datetime, timedelta

ee.Authenticate()

ee.Initialize(project='remotesensing-musa6500')

bbox = ee.Geometry.Polygon([
    [
        [-85.9, 8.0],
        [-85.9, 11.2],
        [-82.5, 11.2],
        [-82.5, 8.0],
        [-85.9, 8.0]
    ]
])

start_date = '2021-07-22'
end_date = '2021-07-28'
start_date = datetime.strptime(start_date, "%Y-%m-%d")
end_date = datetime.strptime(end_date, "%Y-%m-%d")

# Load the datasets

# with the FABDEM layer, the output probability raster 1) is missing tiles and
# 2) shows up semi-transparent in the geemap map for no apparent reason
dem = ee.ImageCollection("projects/sat-io/open-datasets/GLO-30").mosaic().clip(bbox)

# with the hydrosheds DEM, however, this is not the issue
# (uncomment this line to replace the FABDEM layer)
# dem = ee.Image('WWF/HydroSHEDS/03VFDEM').clip(bbox)

slope = ee.Terrain.slope(dem)
landcover = ee.Image("ESA/WorldCover/v100/2020").select('Map').clip(bbox)

stream_dist_proximity_collection = ee.ImageCollection("projects/sat-io/open-datasets/HYDROGRAPHY90/stream-outlet-distance/stream_dist_proximity")\
    .filterBounds(bbox)\
    .mosaic()
stream_dist_proximity = stream_dist_proximity_collection.clip(bbox).rename('stream_distance')

hydro_proj = stream_dist_proximity.projection()

## set time frame
before_start= '2023-09-25'
before_end='2023-10-05'

after_start='2023-10-05'
after_end='2023-10-15'

# SET SAR PARAMETERS (can be left default)

# Polarization (choose either "VH" or "VV")
polarization = "VH"  # or "VV"

# Pass direction (choose either "DESCENDING" or "ASCENDING")
pass_direction = "DESCENDING"  # or "ASCENDING"

# Difference threshold to be applied on the difference image (after flood - before flood)
# It has been chosen by trial and error. Adjust as needed.
difference_threshold = 1.25

# Relative orbit (optional, if you know the relative orbit for your study area)
# relative_orbit = 79

# Rename the selected geometry feature
aoi = bbox

# Load and filter Sentinel-1 GRD data by predefined parameters
collection = ee.ImageCollection('COPERNICUS/S1_GRD') \
    .filter(ee.Filter.eq('instrumentMode', 'IW')) \
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', polarization)) \
    .filter(ee.Filter.eq('orbitProperties_pass', pass_direction)) \
    .filter(ee.Filter.eq('resolution_meters', 10)) \
    .filterBounds(aoi) \
    .select(polarization)

# Select images by predefined dates
before_collection = collection.filterDate(before_start, before_end)
after_collection = collection.filterDate(after_start, after_end)

# Create a mosaic of selected tiles and clip to the study area
before = before_collection.mosaic().clip(aoi)
after = after_collection.mosaic().clip(aoi)

# Apply radar speckle reduction by smoothing
smoothing_radius = 50
before_filtered = before.focal_mean(smoothing_radius, 'circle', 'meters')
after_filtered = after.focal_mean(smoothing_radius, 'circle', 'meters')

# Calculate the difference between the before and after images
difference = after_filtered.divide(before_filtered)

# Apply the predefined difference-threshold and create the flood extent mask
threshold = difference_threshold
difference_binary = difference.gt(threshold)

# Refine the flood result using additional datasets
swater = ee.Image('JRC/GSW1_0/GlobalSurfaceWater').select('seasonality')
swater_mask = swater.gte(10).updateMask(swater.gte(10))
flooded_mask = difference_binary.where(swater_mask, 0)
flooded = flooded_mask.updateMask(flooded_mask)
connections = flooded.connectedPixelCount()
flooded = flooded.updateMask(connections.gte(8))

# Mask out areas with more than 5 percent slope using a Digital Elevation Model
terrain = ee.Algorithms.Terrain(dem)
flooded = flooded.updateMask(slope.lt(5))

# Set the default projection from the hydrography dataset
flooded = flooded.setDefaultProjection(hydro_proj)

# Now, reduce the resolution
flooded_mode = flooded.reduceResolution(
    reducer=ee.Reducer.mode(),
    maxPixels=10000
).reproject(
    crs=hydro_proj
)

# Create a full-area mask, initially marking everything as non-flooded (value 0)
full_area_mask = ee.Image.constant(0).clip(aoi)

# Update the mask to mark flooded areas (value 1)
# Assuming flooded_mode is a binary image with 1 for flooded areas and 0 elsewhere
flood_labeled_image = full_area_mask.where(flooded, 1)

# Now flood_labeled_image contains 1 for flooded areas and 0 for non-flooded areas

combined = (dem.rename("elevation")
        .addBands(landcover.select('Map').rename("landcover"))
        .addBands(slope)
        .addBands(flood_labeled_image.rename("flooded_mask"))
        )

# Get all band names from the combined image
allBandNames = combined.bandNames()

# Remove the class band name ('flooded_full_mask') to get input properties
inputProperties = allBandNames.filter(ee.Filter.neq('item', 'flooded_mask'))

# Perform stratified sampling
stratifiedSample = combined.stratifiedSample(
    numPoints=500,  # Total number of points
    classBand='flooded_mask',  # Band to stratify by
    region=bbox,
    scale=90,
    seed=0
).randomColumn()

# Split into training and testing
training = stratifiedSample.filter(ee.Filter.lt('random', 0.7))
testing = stratifiedSample.filter(ee.Filter.gte('random', 0.7))

# Set up the Random Forest classifier for flood prediction
classifier = ee.Classifier.smileRandomForest(10).train(
    features=training,
    classProperty='flooded_mask',  # Use 'flooded_full_mask' as the class property
    inputProperties=inputProperties  # Dynamically generated input properties
)

# Classify the image
classified = combined.select(inputProperties).classify(classifier)

# Assess accuracy
testAccuracy = testing.classify(classifier).errorMatrix('flooded_mask', 'classification')

# Calculate accuracy
accuracy = testAccuracy.accuracy().getInfo()

# Convert the confusion matrix to an array
confusionMatrixArray = testAccuracy.array().getInfo()

# Calculate recall for the positive class (assuming '1' represents the positive class for flooding)
true_positives = confusionMatrixArray[1][1]  # True positives
false_negatives = confusionMatrixArray[1][0]  # False negatives
false_positives = confusionMatrixArray[0][1]  # False positives (non-flooded incorrectly identified as flooded)
true_negatives = confusionMatrixArray[0][0]  # True negatives (non-flooded correctly identified as non-flooded)
recall = true_positives / (true_positives + false_negatives)
false_positive_rate = false_positives / (false_positives + true_negatives)

print('Confusion Matrix:', confusionMatrixArray)
print('Accuracy:', accuracy)
print('Recall:', recall)
print('False Positive Rate:', false_positive_rate)

Confusion Matrix: [[84, 55], [16, 130]]
Accuracy: 0.7508771929824561
Recall: 0.8904109589041096
False Positive Rate: 0.39568345323741005

# Set up the Random Forest classifier for flood prediction with probability output
classifier = ee.Classifier.smileRandomForest(10).setOutputMode('PROBABILITY').train(
        features=training,
        classProperty='flooded_mask',
        inputProperties=inputProperties
    )

# Classify the image to get probabilities
probabilityImage = combined.select(inputProperties).classify(classifier)

# Visualization parameters for probability with white at the midpoint
vizParamsProbability = {
    'min': 0.,
    'max': 1,
    'palette': colorcet.bmw
}

m = geemap.Map()
m.add("basemap_selector")
m.add("layer_manager")

# Center the map on San Jose, Costa Rica
m.setCenter(-84.0833, 9.9333, 10)
m.addLayer(landcover, {}, 'ESA WorldCover 2020')
m.addLayer(probabilityImage, vizParamsProbability, 'Flood Probability')
m.addLayer(swater_mask, {'palette': 'black'}, 'Permanent Surface Water')
# Display the map
m