6. Analyzing Geospatial Data

6.1. Introduction

6.2. Technical requirements

conda create -n gee python
conda activate gee
conda install -c conda-forge mamba
mamba install -c conda-forge pygis

jupyter lab

# %pip install pygis

import ee
import geemap

geemap.ee_initialize()

6.3. Earth Engine data reductions

6.3.1. List reductions

values = ee.List.sequence(1, 10)
print(values.getInfo())

count = values.reduce(ee.Reducer.count())
print(count.getInfo())  # 10

min_value = values.reduce(ee.Reducer.min())
print(min_value.getInfo())  # 1

max_value = values.reduce(ee.Reducer.max())
print(max_value.getInfo())  # 10

min_max_value = values.reduce(ee.Reducer.minMax())
print(min_max_value.getInfo())

mean_value = values.reduce(ee.Reducer.mean())
print(mean_value.getInfo())  # 5.5

median_value = values.reduce(ee.Reducer.median())
print(median_value.getInfo())  # 5.5

sum_value = values.reduce(ee.Reducer.sum())
print(sum_value.getInfo())  # 55

std_value = values.reduce(ee.Reducer.stdDev())
print(std_value.getInfo())  # 2.8723

6.3.2. ImageCollection reductions

Map = geemap.Map()

# Load an image collection, filtered so it's not too much data.
collection = (
    ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
    .filterDate('2021-01-01', '2021-12-31')
    .filter(ee.Filter.eq('WRS_PATH', 44))
    .filter(ee.Filter.eq('WRS_ROW', 34))
)

# Compute the median in each band, each pixel.
# Band names are B1_median, B2_median, etc.
median = collection.reduce(ee.Reducer.median())

# The output is an Image.  Add it to the map.
vis_param = {'bands': ['B5_median', 'B4_median', 'B3_median'], 'gamma': 2}
Map.setCenter(-122.3355, 37.7924, 8)
Map.addLayer(median, vis_param)
Map

median = collection.median()
print(median.bandNames().getInfo())

6.3.3. Image reductions

Map = geemap.Map()
image = ee.Image('LANDSAT/LC08/C01/T1/LC08_044034_20140318').select(['B4', 'B3', 'B2'])
maxValue = image.reduce(ee.Reducer.max())
Map.centerObject(image, 8)
Map.addLayer(image, {}, 'Original image')
Map.addLayer(maxValue, {'max': 13000}, 'Maximum value image')
Map

6.3.4. FeatureCollection reductions

Map = geemap.Map()
census = ee.FeatureCollection('TIGER/2010/Blocks')
benton = census.filter(
    ee.Filter.And(ee.Filter.eq('statefp10', '41'), ee.Filter.eq('countyfp10', '003'))
)
Map.setCenter(-123.27, 44.57, 13)
Map.addLayer(benton)
Map

# Compute sums of the specified properties.
properties = ['pop10', 'housing10']
sums = benton.filter(ee.Filter.notNull(properties)).reduceColumns(
    **{'reducer': ee.Reducer.sum().repeat(2), 'selectors': properties}
)
sums

print(benton.aggregate_sum('pop10'))  # 85579
print(benton.aggregate_sum('housing10'))  # 36245

benton.aggregate_stats('pop10')

6.4. Image descriptive statistics

Map = geemap.Map()

centroid = ee.Geometry.Point([-122.4439, 37.7538])
image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterBounds(centroid).first()
vis = {'min': 0, 'max': 3000, 'bands': ['B5', 'B4', 'B3']}

Map.centerObject(centroid, 8)
Map.addLayer(image, vis, "Landsat-8")
Map

image.propertyNames()

image.get('CLOUD_COVER')  # 0.05

props = geemap.image_props(image)
props

stats = geemap.image_stats(image, scale=30)
stats

6.5. Zonal statistics with Earth Engine

6.5.1. Zonal statistics

Map = geemap.Map(center=[40, -100], zoom=4)

# Add NASA SRTM
dem = ee.Image('USGS/SRTMGL1_003')
dem_vis = {
    'min': 0,
    'max': 4000,
    'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5'],
}
Map.addLayer(dem, dem_vis, 'SRTM DEM')

# Add 5-year Landsat TOA composite
landsat = ee.Image('LANDSAT/LE7_TOA_5YEAR/1999_2003')
landsat_vis = {'bands': ['B4', 'B3', 'B2'], 'gamma': 1.4}
Map.addLayer(landsat, landsat_vis, "Landsat", False)

# Add US Census States
states = ee.FeatureCollection("TIGER/2018/States")
style = {'fillColor': '00000000'}
Map.addLayer(states.style(**style), {}, 'US States')
Map

out_dem_stats = 'dem_stats.csv'
geemap.zonal_stats(
    dem, states, out_dem_stats, statistics_type='MEAN', scale=1000, return_fc=False
)

out_landsat_stats = 'landsat_stats.csv'
geemap.zonal_stats(
    landsat,
    states,
    out_landsat_stats,
    statistics_type='MEAN',
    scale=1000,
    return_fc=False,
)

6.5.2. Zonal statistics by group

Map = geemap.Map(center=[40, -100], zoom=4)

# Add NLCD data
dataset = ee.Image('USGS/NLCD_RELEASES/2019_REL/NLCD/2019')
landcover = dataset.select('landcover')
Map.addLayer(landcover, {}, 'NLCD 2019')

# Add US census states
states = ee.FeatureCollection("TIGER/2018/States")
style = {'fillColor': '00000000'}
Map.addLayer(states.style(**style), {}, 'US States')

# Add NLCD legend
Map.add_legend(title='NLCD Land Cover', builtin_legend='NLCD')
Map

nlcd_stats = 'nlcd_stats.csv'

geemap.zonal_stats_by_group(
    landcover,
    states,
    nlcd_stats,
    statistics_type='SUM',
    denominator=1e6,
    decimal_places=2,
)

nlcd_stats = 'nlcd_stats_pct.csv'

geemap.zonal_stats_by_group(
    landcover,
    states,
    nlcd_stats,
    statistics_type='PERCENTAGE',
    denominator=1e6,
    decimal_places=2,
)

6.5.3. Zonal statistics with two images

Map = geemap.Map(center=[40, -100], zoom=4)
dem = ee.Image('USGS/3DEP/10m')
vis = {'min': 0, 'max': 4000, 'palette': 'terrain'}
Map.addLayer(dem, vis, 'DEM')
Map

landcover = ee.Image("USGS/NLCD_RELEASES/2019_REL/NLCD/2019").select('landcover')
Map.addLayer(landcover, {}, 'NLCD 2019')
Map.add_legend(title='NLCD Land Cover Classification', builtin_legend='NLCD')

stats = geemap.image_stats_by_zone(dem, landcover, reducer='MEAN')
stats

stats.to_csv('mean.csv', index=False)

geemap.image_stats_by_zone(dem, landcover, out_csv="std.csv", reducer='STD')

6.6. Coordinate grids and fishnets

6.6.1. Creating coordinate grids

lat_grid = geemap.latitude_grid(step=5.0, west=-180, east=180, south=-85, north=85)

Map = geemap.Map()
style = {'fillColor': '00000000'}
Map.addLayer(lat_grid.style(**style), {}, 'Latitude Grid')
Map

df = geemap.ee_to_df(lat_grid)
df

lon_grid = geemap.longitude_grid(step=5.0, west=-180, east=180, south=-85, north=85)

Map = geemap.Map()
style = {'fillColor': '00000000'}
Map.addLayer(lon_grid.style(**style), {}, 'Longitude Grid')
Map

grid = geemap.latlon_grid(
    lat_step=10, lon_step=10, west=-180, east=180, south=-85, north=85
)

Map = geemap.Map()
style = {'fillColor': '00000000'}
Map.addLayer(grid.style(**style), {}, 'Coordinate Grid')
Map

6.6.2. Creating fishnets

Map = geemap.Map()
Map

roi = Map.user_roi

if roi is None:
    roi = ee.Geometry.BBox(-112.8089, 33.7306, -88.5951, 46.6244)
    Map.addLayer(roi, {}, 'ROI')
    Map.user_roi = None

Map.centerObject(roi)

fishnet = geemap.fishnet(roi, h_interval=2.0, v_interval=2.0, delta=1)
style = {'color': 'blue', 'fillColor': '00000000'}
Map.addLayer(fishnet.style(**style), {}, 'Fishnet')

Map = geemap.Map()
Map

roi = Map.user_roi

if roi is None:
    roi = ee.Geometry.Polygon(
        [
            [
                [-64.602356, -1.127399],
                [-68.821106, -12.625598],
                [-60.647278, -22.498601],
                [-47.815247, -21.111406],
                [-43.860168, -8.913564],
                [-54.582825, -0.775886],
                [-60.823059, 0.454555],
                [-64.602356, -1.127399],
            ]
        ]
    )
    Map.addLayer(roi, {}, 'ROI')

Map.centerObject(roi)
Map

fishnet = geemap.fishnet(roi, rows=6, cols=8, delta=1)
style = {'color': 'blue', 'fillColor': '00000000'}
Map.addLayer(fishnet.style(**style), {}, 'Fishnet')

6.7. Extracting pixel values

6.7.1. Extracting values to points

Map = geemap.Map(center=[40, -100], zoom=4)

dem = ee.Image('USGS/SRTMGL1_003')
landsat7 = ee.Image('LANDSAT/LE7_TOA_5YEAR/1999_2003')

vis_params = {
    'min': 0,
    'max': 4000,
    'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5'],
}

Map.addLayer(
    landsat7,
    {'bands': ['B4', 'B3', 'B2'], 'min': 20, 'max': 200, 'gamma': 2},
    'Landsat 7',
)
Map.addLayer(dem, vis_params, 'SRTM DEM', True, 1)
Map

in_shp = 'us_cities.shp'
url = 'https://github.com/giswqs/data/raw/main/us/us_cities.zip'
geemap.download_file(url)

in_fc = geemap.shp_to_ee(in_shp)
Map.addLayer(in_fc, {}, 'Cities')

geemap.extract_values_to_points(in_fc, dem, out_fc="dem.shp")

geemap.shp_to_gdf("dem.shp")

geemap.extract_values_to_points(in_fc, landsat7, 'landsat.csv')

geemap.csv_to_df('landsat.csv')

6.7.2. Extracting pixel values along a transect

Map = geemap.Map(center=[40, -100], zoom=4)
Map.add_basemap("TERRAIN")

image = ee.Image('USGS/SRTMGL1_003')
vis_params = {
    'min': 0,
    'max': 4000,
    'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5'],
}
Map.addLayer(image, vis_params, 'SRTM DEM', True, 0.5)
Map

line = Map.user_roi
if line is None:
    line = ee.Geometry.LineString(
        [[-120.2232, 36.3148], [-118.9269, 36.7121], [-117.2022, 36.7562]]
    )
    Map.addLayer(line, {}, "ROI")
Map.centerObject(line)

reducer = 'mean'
transect = geemap.extract_transect(
    image, line, n_segments=100, reducer=reducer, to_pandas=True
)
transect

geemap.line_chart(
    data=transect,
    x='distance',
    y='mean',
    markers=True,
    x_label='Distance (m)',
    y_label='Elevation (m)',
    height=400,
)

transect.to_csv('transect.csv')

6.7.3. Interactive region reduction

Map = geemap.Map()

collection = (
    ee.ImageCollection('MODIS/061/MOD13A2')
    .filterDate('2015-01-01', '2019-12-31')
    .select('NDVI')
)

image = collection.toBands()

ndvi_vis = {
    'min': 0.0,
    'max': 9000.0,
    'palette': 'ndvi',
}

Map.addLayer(image, {}, 'MODIS NDVI Time-series')
Map.addLayer(image.select(0), ndvi_vis, 'First image')

Map

dates = geemap.image_dates(collection).getInfo()
dates

len(dates)

Map.set_plot_options(add_marker_cluster=True)
Map.roi_reducer = ee.Reducer.mean()
Map

Map.extract_values_to_points('ndvi.csv')

6.8. Mapping available image count

collection = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
image = geemap.image_count(
    collection, region=None, start_date='2021-01-01', end_date='2022-01-01', clip=False
)

Map = geemap.Map()
vis = {'min': 0, 'max': 60, 'palette': 'coolwarm'}
Map.addLayer(image, vis, 'Image Count')
Map.add_colorbar(vis, label='Landsat 8 Image Count')

countries = ee.FeatureCollection(geemap.examples.get_ee_path('countries'))
style = {"color": "00000088", "width": 1, "fillColor": "00000000"}
Map.addLayer(countries.style(**style), {}, "Countries")
Map

6.9. Cloud-free composites

Map = geemap.Map()

collection = ee.ImageCollection('LANDSAT/LC08/C02/T1').filterDate(
    '2021-01-01', '2022-01-01'
)

composite = ee.Algorithms.Landsat.simpleComposite(collection)

vis_params = {'bands': ['B5', 'B4', 'B3'], 'max': 128}

Map.setCenter(-122.3578, 37.7726, 10)
Map.addLayer(composite, vis_params, 'TOA composite')
Map

customComposite = ee.Algorithms.Landsat.simpleComposite(
    **{'collection': collection, 'percentile': 30, 'cloudScoreRange': 5}
)

Map.addLayer(customComposite, vis_params, 'Custom TOA composite')
Map.setCenter(-105.4317, 52.5536, 11)

vis_params = [
    {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 128},
    {'bands': ['B5', 'B4', 'B3'], 'min': 0, 'max': 128},
    {'bands': ['B7', 'B6', 'B4'], 'min': 0, 'max': 128},
    {'bands': ['B6', 'B5', 'B2'], 'min': 0, 'max': 128},
]

labels = [
    'Natural Color (4, 3, 2)',
    'Color Infrared (5, 4, 3)',
    'Short-Wave Infrared (7, 6 4)',
    'Agriculture (6, 5, 2)',
]

geemap.linked_maps(
    rows=2,
    cols=2,
    height="300px",
    center=[37.7726, -122.1578],
    zoom=9,
    ee_objects=[composite],
    vis_params=vis_params,
    labels=labels,
    label_position="topright",
)

6.10. Quality mosaicking

Map = geemap.Map(center=[40, -100], zoom=4)
countries = ee.FeatureCollection(geemap.examples.get_ee_path('countries'))
roi = countries.filter(ee.Filter.eq('ISO_A3', 'USA'))
Map.addLayer(roi, {}, 'roi')
Map

start_date = '2020-01-01'
end_date = '2021-01-01'
collection = (
    ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
    .filterBounds(roi)
    .filterDate(start_date, end_date)
)

median = collection.median()
vis_rgb = {
    'bands': ['B4', 'B3', 'B2'],
    'min': 0,
    'max': 0.4,
}
Map.addLayer(median, vis_rgb, 'Median')
Map

def add_ndvi(image):
    ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI')
    return image.addBands(ndvi)

def add_time(image):
    date = ee.Date(image.date())

    img_date = ee.Number.parse(date.format('YYYYMMdd'))
    image = image.addBands(ee.Image(img_date).rename('date').toInt())

    img_month = ee.Number.parse(date.format('M'))
    image = image.addBands(ee.Image(img_month).rename('month').toInt())

    img_doy = ee.Number.parse(date.format('D'))
    image = image.addBands(ee.Image(img_doy).rename('doy').toInt())

    return image

images = collection.map(add_ndvi).map(add_time)

greenest = images.qualityMosaic('NDVI')

greenest.bandNames()

ndvi = greenest.select('NDVI')
vis_ndvi = {'min': 0, 'max': 1, 'palette': 'ndvi'}
Map.addLayer(ndvi, vis_ndvi, 'NDVI')
Map.add_colorbar(vis_ndvi, label='NDVI', layer_name='NDVI')
Map

Map.addLayer(greenest, vis_rgb, 'Greenest pixel')

vis_month = {'palette': ['red', 'blue'], 'min': 1, 'max': 12}
Map.addLayer(greenest.select('month'), vis_month, 'Greenest month')
Map.add_colorbar(vis_month, label='Month', layer_name='Greenest month')

vis_doy = {'palette': ['brown', 'green'], 'min': 1, 'max': 365}
Map.addLayer(greenest.select('doy'), vis_doy, 'Greenest doy')
Map.add_colorbar(vis_doy, label='Day of year', layer_name='Greenest doy')

6.11. Interactive charts

6.11.1. Chart Overview

6.11.2. Data table charts

data = geemap.examples.get_path('countries.geojson')
df = geemap.geojson_to_df(data)
df.head()

geemap.bar_chart(
    data=df,
    x='NAME',
    y='POP_EST',
    x_label='Country',
    y_label='Population',
    descending=True,
    max_rows=30,
    title='World Population',
    height=500,
    layout_args={'title_x': 0.5, 'title_y': 0.85},
)

geemap.pie_chart(
    data=df,
    names='NAME',
    values='POP_EST',
    max_rows=30,
    height=600,
    title='World Population',
    legend_title='Country',
    layout_args={'title_x': 0.47, 'title_y': 0.87},
)

data = geemap.examples.get_path('life_exp.csv')
df = geemap.csv_to_df(data)
df = df[df['continent'] == 'Oceania']
df.head()

geemap.line_chart(
    df,
    x='year',
    y='lifeExp',
    color='country',
    x_label='Year',
    y_label='Life expectancy',
    legend_title='Country',
    height=400,
    markers=True,
)

6.11.3. Earth Engine object charts

import geemap.chart as chart

Map = geemap.Map(center=[40, -100], zoom=4)
collection = ee.FeatureCollection('projects/google/charts_feature_example')
Map.addLayer(collection, {}, "Ecoregions")
Map

6.11.3.1. Chart by feature

features = collection.select('[0-9][0-9]_tmean|label')
df = geemap.ee_to_df(features, sort_columns=True)
df

xProperty = "label"
yProperties = df.columns[:12]

labels = [
    'Jan',
    'Feb',
    'Mar',
    'Apr',
    'May',
    'Jun',
    'Jul',
    'Aug',
    'Sep',
    'Oct',
    'Nov',
    'Dec',
]
colors = [
    '#604791',
    '#1d6b99',
    '#39a8a7',
    '#0f8755',
    '#76b349',
    '#f0af07',
    '#e37d05',
    '#cf513e',
    '#96356f',
    '#724173',
    '#9c4f97',
    '#696969',
]
title = "Average Monthly Temperature by Ecoregion"
xlabel = "Ecoregion"
ylabel = "Temperature"

options = {
    "labels": labels,
    "colors": colors,
    "title": title,
    "xlabel": xlabel,
    "ylabel": ylabel,
    "legend_location": "top-left",
    "height": "500px",
}

chart.feature_byFeature(features, xProperty, yProperties, **options)

6.11.3.2. Chart by property

features = collection.select('[0-9][0-9]_ppt|label')
df = geemap.ee_to_df(features, sort_columns=True)
df

keys = df.columns[:12]
values = [
    'Jan',
    'Feb',
    'Mar',
    'Apr',
    'May',
    'Jun',
    'Jul',
    'Aug',
    'Sep',
    'Oct',
    'Nov',
    'Dec',
]
xProperties = dict(zip(keys, values))
seriesProperty = "label"

options = {
    'title': "Average Ecoregion Precipitation by Month",
    'colors': ['#f0af07', '#0f8755', '#76b349'],
    'xlabel': "Month",
    'ylabel': "Precipitation (mm)",
    'legend_location': "top-left",
    "height": "500px",
}

chart.feature_byProperty(features, xProperties, seriesProperty, **options)

6.11.3.3. Feature histograms

source = ee.ImageCollection('OREGONSTATE/PRISM/Norm81m').toBands()
region = ee.Geometry.Rectangle(-123.41, 40.43, -116.38, 45.14)
samples = source.sample(region, 5000)
prop = '07_ppt'

options = {
    "title": 'July Precipitation Distribution for NW USA',
    "xlabel": 'Precipitation (mm)',
    "ylabel": 'Pixel count',
    "colors": ['#1d6b99'],
}

chart.feature_histogram(samples, prop, **options)

chart.feature_histogram(samples, prop, maxBuckets=30, **options)

chart.feature_histogram(samples, prop, minBucketWidth=0.5, **options)

chart.feature_histogram(samples, prop, minBucketWidth=3, maxBuckets=30, **options)

6.12. Unsupervised classification

Map = geemap.Map()

point = ee.Geometry.Point([-88.0664, 41.9411])

image = (
    ee.ImageCollection('LANDSAT/LC09/C02/T1_L2')
    .filterBounds(point)
    .filterDate('2022-01-01', '2022-12-31')
    .sort('CLOUD_COVER')
    .first()
    .select('SR_B[1-7]')
)

region = image.geometry()
image = image.multiply(0.0000275).add(-0.2).set(image.toDictionary())
vis_params = {'min': 0, 'max': 0.3, 'bands': ['SR_B5', 'SR_B4', 'SR_B3']}

Map.centerObject(region, 8)
Map.addLayer(image, vis_params, "Landsat-9")
Map

geemap.get_info(image)

image.get('DATE_ACQUIRED').getInfo()

image.get('CLOUD_COVER').getInfo()

training = image.sample(
    **{
        # "region": region,
        'scale': 30,
        'numPixels': 5000,
        'seed': 0,
        'geometries': True,  # Set this to False to ignore geometries
    }
)

Map.addLayer(training, {}, 'Training samples')
Map

geemap.ee_to_df(training.limit(5))

n_clusters = 5
clusterer = ee.Clusterer.wekaKMeans(n_clusters).train(training)

result = image.cluster(clusterer)
Map.addLayer(result.randomVisualizer(), {}, 'clusters')
Map

legend_dict = {
    'Open Water': '#466b9f',
    'Developed, High Intensity': '#ab0000',
    'Developed, Low Intensity': '#d99282',
    'Forest': '#1c5f2c',
    'Cropland': '#ab6c28'

}

palette = list(legend_dict.values())

Map.addLayer(
    result, {'min': 0, 'max': 4, 'palette': palette}, 'Labelled clusters'
)
Map.add_legend(title='Land Cover Type',legend_dict=legend_dict , position='bottomright')
Map

geemap.download_ee_image(image, filename='unsupervised.tif', region=region, scale=90)

6.13. Supervised classification

Map = geemap.Map()
point = ee.Geometry.Point([-122.4439, 37.7538])

image = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(point)
    .filterDate('2019-01-01', '2020-01-01')
    .sort('CLOUD_COVER')
    .first()
    .select('SR_B[1-7]')
)

image = image.multiply(0.0000275).add(-0.2).set(image.toDictionary())
vis_params = {'min': 0, 'max': 0.3, 'bands': ['SR_B5', 'SR_B4', 'SR_B3']}

Map.centerObject(point, 8)
Map.addLayer(image, vis_params, "Landsat-8")
Map

geemap.get_info(image)

image.get('DATE_ACQUIRED').getInfo()

image.get('CLOUD_COVER').getInfo()

nlcd = ee.Image('USGS/NLCD_RELEASES/2019_REL/NLCD/2019')
landcover = nlcd.select('landcover').clip(image.geometry())
Map.addLayer(landcover, {}, 'NLCD Landcover')
Map

points = landcover.sample(
    **{
        'region': image.geometry(),
        'scale': 30,
        'numPixels': 5000,
        'seed': 0,
        'geometries': True,
    }
)

Map.addLayer(points, {}, 'training', False)

print(points.size().getInfo())

bands = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7']
label = 'landcover'
features = image.select(bands).sampleRegions(
    **{'collection': points, 'properties': [label], 'scale': 30}
)

geemap.ee_to_df(features.limit(5))

params = {

    'features': features,
    'classProperty': label,
    'inputProperties': bands,

}
classifier = ee.Classifier.smileCart(maxNodes=None).train(**params)

classified = image.select(bands).classify(classifier).rename('landcover')
Map.addLayer(classified.randomVisualizer(), {}, 'Classified')
Map

geemap.get_info(nlcd)

class_values = nlcd.get('landcover_class_values')
class_palette = nlcd.get('landcover_class_palette')
classified = classified.set({
    'landcover_class_values': class_values,
    'landcover_class_palette': class_palette
})

Map.addLayer(classified, {}, 'Land cover')
Map.add_legend(title="Land cover type", builtin_legend='NLCD')
Map

geemap.download_ee_image(
    landcover,
    filename='supervised.tif',
    region=image.geometry(),
    scale=30
    )

6.14. Accuracy assessment

Map = geemap.Map()
point = ee.Geometry.Point([-122.4439, 37.7538])

img = (
    ee.ImageCollection('COPERNICUS/S2_SR')
    .filterBounds(point)
    .filterDate('2020-01-01', '2021-01-01')
    .sort('CLOUDY_PIXEL_PERCENTAGE')
    .first()
    .select('B.*')
)

vis_params = {'min': 100, 'max': 3500, 'bands': ['B11',  'B8',  'B3']}

Map.centerObject(point, 9)
Map.addLayer(img, vis_params, "Sentinel-2")
Map

lc = ee.Image('ESA/WorldCover/v100/2020')
classValues = [10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100]
remapValues = ee.List.sequence(0, 10)
label = 'lc'
lc = lc.remap(classValues, remapValues).rename(label).toByte()

sample = img.addBands(lc).stratifiedSample(**{
  'numPoints': 100,
  'classBand': label,
  'region': img.geometry(),
  'scale': 10,
  'geometries': True
})

sample = sample.randomColumn()
trainingSample = sample.filter('random <= 0.8')
validationSample = sample.filter('random > 0.8')

trainedClassifier = ee.Classifier.smileRandomForest(numberOfTrees=10).train(**{
  'features': trainingSample,
  'classProperty': label,
  'inputProperties': img.bandNames()
})

print('Results of trained classifier', trainedClassifier.explain().getInfo())

trainAccuracy = trainedClassifier.confusionMatrix()
trainAccuracy.getInfo()

trainAccuracy.accuracy().getInfo()

trainAccuracy.kappa().getInfo()

validationSample = validationSample.classify(trainedClassifier)
validationAccuracy = validationSample.errorMatrix(label, 'classification')
validationAccuracy.getInfo()

validationAccuracy.accuracy().getInfo()

validationAccuracy.producersAccuracy().getInfo()

validationAccuracy.consumersAccuracy().getInfo()

import csv

with open("training.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerows(trainAccuracy.getInfo())

with open("validation.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerows(validationAccuracy.getInfo())

imgClassified = img.classify(trainedClassifier)

classVis = {
  'min': 0,
  'max': 10,
  'palette': ['006400' ,'ffbb22', 'ffff4c', 'f096ff', 'fa0000', 'b4b4b4',
            'f0f0f0', '0064c8', '0096a0', '00cf75', 'fae6a0']
}
Map.addLayer(lc, classVis, 'ESA Land Cover', False)
Map.addLayer(imgClassified, classVis, 'Classified')
Map.addLayer(trainingSample, {'color': 'black'}, 'Training sample')
Map.addLayer(validationSample, {'color': 'white'}, 'Validation sample')
Map.add_legend(title='Land Cover Type', builtin_legend='ESA_WorldCover')
Map.centerObject(img)
Map

6.15. Using locally trained machine learning models

import pandas as pd
from geemap import ml
from sklearn import ensemble

6.15.1. Train a model locally using scikit-learn

url = "https://raw.githubusercontent.com/gee-community/geemap/master/examples/data/rf_example.csv"
df = pd.read_csv(url)
df

feature_names = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']
label = "landcover"

X = df[feature_names]
y = df[label]
n_trees = 10
rf = ensemble.RandomForestClassifier(n_trees).fit(X, y)

6.15.2. Convert a sklearn classifier object to a list of strings

trees = ml.rf_to_strings(rf, feature_names)

print(len(trees))

print(trees[0])

6.15.3. Convert sklearn classifier to GEE classifier

ee_classifier = ml.strings_to_classifier(trees)
ee_classifier.getInfo()

6.15.4. Classify image using GEE classifier

# Make a cloud-free Landsat 8 TOA composite (from raw imagery).
l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1')

image = ee.Algorithms.Landsat.simpleComposite(
    collection=l8.filterDate('2018-01-01', '2018-12-31'), asFloat=True
)

classified = image.select(feature_names).classify(ee_classifier)

Map = geemap.Map(center=(37.75, -122.25), zoom=11)

Map.addLayer(
    image,
    {"bands": ['B7', 'B5', 'B3'], "min": 0.05, "max": 0.55, "gamma": 1.5},
    'image',
)
Map.addLayer(
    classified,
    {"min": 0, "max": 2, "palette": ['red', 'green', 'blue']},
    'classification',
)
Map

6.15.5. Save trees to the cloud

user_id = geemap.ee_user_id()
asset_id = user_id + "/random_forest_strings_test"
asset_id

ml.export_trees_to_fc(trees, asset_id)

rf_fc = ee.FeatureCollection(asset_id)
another_classifier = ml.fc_to_classifier(rf_fc)
classified = image.select(feature_names).classify(another_classifier)

6.15.6. Save trees locally

out_csv = "trees.csv"
ml.trees_to_csv(trees, out_csv)
another_classifier = ml.csv_to_classifier(out_csv)
classified = image.select(feature_names).classify(another_classifier)

6.16. Sankey diagrams

import sankee

sankee.datasets.LCMS_LC.sankify(
    years=[1990, 2000, 2010, 2020],
    region=ee.Geometry.Point([-122.192688, 46.25917]).buffer(2000),
    max_classes=3,
    title="Mount St. Helens Recovery",
)

Map = geemap.Map(height=650)
Map

6.17. Summary