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 geemap pygis
jupyter lab

Open in Colab

# 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()

# Load an image and select some bands of interest.
image = ee.Image('LANDSAT/LC08/C01/T1/LC08_044034_20140318').select(['B4', 'B3', 'B2'])

# Reduce the image to get a one-band maximum value image.
maxValue = image.reduce(ee.Reducer.max())

# Display the result.
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()

# Load US census data as a FeatureCollection.
census = ee.FeatureCollection('TIGER/2010/Blocks')

# Filter the collection to include only Benton County, OR.
benton = census.filter(
    ee.Filter.And(ee.Filter.eq('statefp10', '41'), ee.Filter.eq('countyfp10', '003'))
)

# Display Benton County census blocks.
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}
)

# Print the resultant Dictionary.
print(sums.getInfo())
print(benton.aggregate_sum('pop10').getInfo())  # 85579
print(benton.aggregate_sum('housing10').getInfo())  # 36245
benton.aggregate_stats('pop10').getInfo()

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().getInfo()
image.get('CLOUD_COVER').getInfo()  # 0.05
props = geemap.image_props(image)
props.getInfo()
stats = geemap.image_stats(image, scale=30)
stats.getInfo()

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 fishets#

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(data, {}, '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')
print(greenest.bandNames().getInfo())
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. Creating training samples#

Map = geemap.Map()
Map
if Map.user_rois is not None:
    training_samples = Map.user_rois
    print(training_samples.getInfo())

6.13. Unsupervised classification#

6.13.1. Unsupervised classification algorithms#

6.13.2. Add data to the map#

Map = geemap.Map()

# point = ee.Geometry.Point([-122.4439, 37.7538])
point = ee.Geometry.Point([-87.7719, 41.8799])

image = (
    ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
    .filterBounds(point)
    .filterDate('2019-01-01', '2019-12-31')
    .sort('CLOUD_COVER')
    .first()
    .select('B[1-7]')
)

vis_params = {'min': 0, 'max': 3000, 'bands': ['B5', 'B4', 'B3']}

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

6.13.3. Check image properties#

props = geemap.image_props(image)
props.getInfo()
props.get('IMAGE_DATE').getInfo()
props.get('CLOUD_COVER').getInfo()

6.13.4. Create training samples#

# region = Map.user_roi
# region = ee.Geometry.Rectangle([-122.6003, 37.4831, -121.8036, 37.8288])
# region = ee.Geometry.Point([-122.4439, 37.7538]).buffer(10000)
# Make the training dataset.
training = image.sample(
    **{
        #     'region': region,
        'scale': 30,
        'numPixels': 5000,
        'seed': 0,
        'geometries': True,  # Set this to False to ignore geometries
    }
)

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

6.13.5. Train the clusterer#

# Instantiate the clusterer and train it.
n_clusters = 5
clusterer = ee.Clusterer.wekaKMeans(n_clusters).train(training)

6.13.6. Classify the image#

# Cluster the input using the trained clusterer.
result = image.cluster(clusterer)

# # Display the clusters with random colors.
Map.addLayer(result.randomVisualizer(), {}, 'clusters')
Map

6.13.7. Label the clusters#

legend_keys = ['One', 'Two', 'Three', 'Four', 'ect']
legend_colors = ['#8DD3C7', '#FFFFB3', '#BEBADA', '#FB8072', '#80B1D3']

# Reclassify the map
result = result.remap([0, 1, 2, 3, 4], [1, 2, 3, 4, 5])

Map.addLayer(
    result, {'min': 1, 'max': 5, 'palette': legend_colors}, 'Labelled clusters'
)
Map.add_legend(
    legend_keys=legend_keys, legend_colors=legend_colors, position='bottomright'
)
Map

6.13.8. Visualize results#

print('Change layer opacity:')
cluster_layer = Map.layers[-1]
cluster_layer.interact(opacity=(0, 1, 0.1))

6.13.9. Export results#

geemap.ee_export_image(result, filename="cluster.tif", scale=90)
geemap.ee_export_image_to_drive(
    result, description='clusters', folder='export', scale=90
)

6.14. Supervised classification#

6.14.1. Supervised classification algorithms#

6.14.2. Add data to the map#

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

image = (
    ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
    .filterBounds(point)
    .filterDate('2016-01-01', '2016-12-31')
    .sort('CLOUD_COVER')
    .first()
    .select('B[1-7]')
)

vis_params = {'min': 0, 'max': 3000, 'bands': ['B5', 'B4', 'B3']}

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

6.14.3. Check image properties#

ee.Date(image.get('system:time_start')).format('YYYY-MM-dd').getInfo()
image.get('CLOUD_COVER').getInfo()

6.14.4. Create training samples#

# region = Map.user_roi
# region = ee.Geometry.Rectangle([-122.6003, 37.4831, -121.8036, 37.8288])
# region = ee.Geometry.Point([-122.4439, 37.7538]).buffer(10000)

nlcd = ee.Image('USGS/NLCD/NLCD2016').select('landcover').clip(image.geometry())
Map.addLayer(nlcd, {}, 'NLCD')
Map
# Make the training dataset.
points = nlcd.sample(
    **{
        'region': image.geometry(),
        'scale': 30,
        'numPixels': 5000,
        'seed': 0,
        'geometries': True,  # Set this to False to ignore geometries
    }
)

Map.addLayer(points, {}, 'training', False)
print(points.size().getInfo())
print(points.first().getInfo())

6.14.5. Train the classifier#

# Use these bands for prediction.
bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7']


# This property of the table stores the land cover labels.
label = 'landcover'

# Overlay the points on the imagery to get training.
training = image.select(bands).sampleRegions(
    **{'collection': points, 'properties': [label], 'scale': 30}
)

# Train a CART classifier with default parameters.
trained = ee.Classifier.smileCart().train(training, label, bands)
print(training.first().getInfo())

6.14.6. Classify the image#

# Classify the image with the same bands used for training.
result = image.select(bands).classify(trained)

# # Display the clusters with random colors.
Map.addLayer(result.randomVisualizer(), {}, 'classified')
Map

6.14.7. Render categorical map#

class_values = nlcd.get('landcover_class_values').getInfo()
class_values
class_palette = nlcd.get('landcover_class_palette').getInfo()
class_palette
landcover = result.set('classification_class_values', class_values)
landcover = landcover.set('classification_class_palette', class_palette)
Map.addLayer(landcover, {}, 'Land cover')
Map

6.14.8. Visualize results#

print('Change layer opacity:')
cluster_layer = Map.layers[-1]
cluster_layer.interact(opacity=(0, 1, 0.1))

6.14.9. Add a legend to the map#

Map.add_legend(builtin_legend='NLCD')
Map

6.14.10. Export results#

geemap.ee_export_image(landcover, filename='landcover.tif', scale=900)
geemap.ee_export_image_to_drive(
    landcover, description='landcover', folder='export', scale=900
)

6.15. Accuracy assessment#

6.15.1. Supervised classification algorithms#

6.15.2. Add data to the map#

Map = geemap.Map()
NLCD2016 = ee.Image('USGS/NLCD/NLCD2016').select('landcover')
Map.addLayer(NLCD2016, {}, 'NLCD 2016')
Map
NLCD_metadata = ee.FeatureCollection("users/giswqs/landcover/NLCD2016_metadata")
Map.addLayer(NLCD_metadata, {}, 'NLCD Metadata')
# point = ee.Geometry.Point([-122.4439, 37.7538])  # Sanfrancisco, CA
# point = ee.Geometry.Point([-83.9293, 36.0526])   # Knoxville, TN
point = ee.Geometry.Point([-88.3070, 41.7471])  # Chicago, IL
metadata = NLCD_metadata.filterBounds(point).first()
region = metadata.geometry()
metadata.get('2016on_bas').getInfo()
doy = metadata.get('2016on_bas').getInfo().replace('LC08_', '')
doy
ee.Date.parse('YYYYDDD', doy).format('YYYY-MM-dd').getInfo()
start_date = ee.Date.parse('YYYYDDD', doy)
end_date = start_date.advance(1, 'day')
image = (
    ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
    .filterBounds(point)
    .filterDate(start_date, end_date)
    .first()
    .select('B[1-7]')
    .clip(region)
)

vis_params = {'min': 0, 'max': 3000, 'bands': ['B5', 'B4', 'B3']}

Map.centerObject(point, 8)
Map.addLayer(image, vis_params, "Landsat-8")
Map
nlcd_raw = NLCD2016.clip(region)
Map.addLayer(nlcd_raw, {}, 'NLCD')

6.15.3. Prepare for consecutive class labels#

raw_class_values = nlcd_raw.get('landcover_class_values').getInfo()
print(raw_class_values)
n_classes = len(raw_class_values)
new_class_values = list(range(0, n_classes))
new_class_values
class_palette = nlcd_raw.get('landcover_class_palette').getInfo()
print(class_palette)
nlcd = nlcd_raw.remap(raw_class_values, new_class_values).select(
    ['remapped'], ['landcover']
)
nlcd = nlcd.set('landcover_class_values', new_class_values)
nlcd = nlcd.set('landcover_class_palette', class_palette)
Map.addLayer(nlcd, {}, 'NLCD')
Map

6.15.4. Make training data#

# Make the training dataset.
points = nlcd.sample(
    **{
        'region': region,
        'scale': 30,
        'numPixels': 5000,
        'seed': 0,
        'geometries': True,  # Set this to False to ignore geometries
    }
)

Map.addLayer(points, {}, 'training', False)
print(points.size().getInfo())
print(points.first().getInfo())

6.15.5. Split training and testing#

# Use these bands for prediction.
bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7']

# This property of the table stores the land cover labels.
label = 'landcover'

# Overlay the points on the imagery to get training.
sample = image.select(bands).sampleRegions(
    **{'collection': points, 'properties': [label], 'scale': 30}
)

# Adds a column of deterministic pseudorandom numbers.
sample = sample.randomColumn()

split = 0.7

training = sample.filter(ee.Filter.lt('random', split))
validation = sample.filter(ee.Filter.gte('random', split))
training.first().getInfo()
validation.first().getInfo()

6.15.6. Train the classifier#

classifier = ee.Classifier.smileRandomForest(10).train(training, label, bands)

6.15.7. Classify the image#

# Classify the image with the same bands used for training.
result = image.select(bands).classify(classifier)

# # Display the clusters with random colors.
Map.addLayer(result.randomVisualizer(), {}, 'classfied')
Map

6.15.8. Render categorical map#

class_values = nlcd.get('landcover_class_values').getInfo()
print(class_values)
class_palette = nlcd.get('landcover_class_palette').getInfo()
print(class_palette)
landcover = result.set('classification_class_values', class_values)
landcover = landcover.set('classification_class_palette', class_palette)
Map.addLayer(landcover, {}, 'Land cover')
Map

6.15.9. Visualize results#

print('Change layer opacity:')
cluster_layer = Map.layers[-1]
cluster_layer.interact(opacity=(0, 1, 0.1))

6.15.10. Add a legend to the map#

Map.add_legend(builtin_legend='NLCD')
Map

6.15.11. Accuracy assessment#

6.15.11.1. Training dataset#

train_accuracy = classifier.confusionMatrix()
train_accuracy.getInfo()
train_accuracy.accuracy().getInfo()
train_accuracy.kappa().getInfo()
train_accuracy.producersAccuracy().getInfo()
train_accuracy.consumersAccuracy().getInfo()
6.15.11.1.1. Validation dataset#
validated = validation.classify(classifier)
validated.first().getInfo()
test_accuracy = validated.errorMatrix('landcover', 'classification')
test_accuracy.getInfo()
test_accuracy.accuracy().getInfo()
test_accuracy.kappa().getInfo()
test_accuracy.producersAccuracy().getInfo()
test_accuracy.consumersAccuracy().getInfo()

6.15.12. Download confusion matrix#

import csv
import os

out_dir = os.getcwd()
training_csv = os.path.join(out_dir, 'train_accuracy.csv')
testing_csv = os.path.join(out_dir, 'test_accuracy.csv')

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

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

6.15.13. Reclassify land cover map#

landcover = landcover.remap(new_class_values, raw_class_values).select(
    ['remapped'], ['classification']
)
landcover = landcover.set('classification_class_values', raw_class_values)
landcover = landcover.set('classification_class_palette', class_palette)
Map.addLayer(landcover, {}, 'Final land cover')
Map

6.15.14. Export results#

import os

out_dir = os.getcwd()
out_file = os.path.join(out_dir, 'landcover.tif')
geemap.ee_export_image(landcover, filename=out_file, scale=900)
geemap.ee_export_image_to_drive(
    landcover, description='landcover', folder='export', scale=900
)

6.16. Using locally trained machine learning models#

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

6.16.1. Train a model locally using scikit-learn#

# read the feature table to train our RandomForest model
# data taken from ee.FeatureCollection('GOOGLE/EE/DEMOS/demo_landcover_labels')

url = "https://raw.githubusercontent.com/giswqs/geemap/master/examples/data/rf_example.csv"
df = pd.read_csv(url)
df
# specify the names of the features (i.e. band names) and label
# feature names used to extract out features and define what bands

feature_names = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']
label = "landcover"
# get the features and labels into separate variables
X = df[feature_names]
y = df[label]
# create a classifier and fit
n_trees = 10
rf = ensemble.RandomForestClassifier(n_trees).fit(X, y)

6.16.2. Convert a sklearn classifier object to a list of strings#

# convert the estimator into a list of strings
# this function also works with the ensemble.ExtraTrees estimator
trees = ml.rf_to_strings(rf, feature_names)
# print the first tree to see the result
print(trees[0])
print(trees[1])
# number of trees we converted should equal the number of trees we defined for the model
len(trees) == n_trees

6.16.3. Convert sklearn classifier to GEE classifier#

# create a ee classifier to use with ee objects from the trees
ee_classifier = ml.strings_to_classifier(trees)
# ee_classifier.getInfo()

6.16.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
)
# classify the image using the classifier we created from the local training
# note: here we select the feature_names from the image that way the classifier knows which bands to use
classified = image.select(feature_names).classify(ee_classifier)
# display results
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.16.5. Save trees to the cloud#

user_id = geemap.ee_user_id()
user_id
# specify asset id where to save trees
# be sure to change <user_name> to your ee user name
asset_id = user_id + "/random_forest_strings_test"
asset_id
# kick off an export process so it will be saved to the ee asset
ml.export_trees_to_fc(trees, asset_id)

# this will kick off an export task, so wait a few minutes before moving on
# read the exported tree feature collection
rf_fc = ee.FeatureCollection(asset_id)

# convert it to a classifier, very similar to the `ml.trees_to_classifier` function
another_classifier = ml.fc_to_classifier(rf_fc)

# classify the image again but with the classifier from the persisted trees
classified = image.select(feature_names).classify(another_classifier)
# display results
# we should get the exact same results as before
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.16.6. Save trees locally#

import os

out_csv = os.path.abspath("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)
# display results
# we should get the exact same results as before
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.17. Global land cover area#

Map = geemap.Map()
dataset = ee.ImageCollection("ESA/WorldCover/v100").first()
Map.addLayer(dataset, {'bands': ['Map']}, 'ESA Land Cover')
Map.add_legend(builtin_legend='ESA_WorldCover')
Map
df = geemap.image_area_by_group(
    dataset, scale=1000, denominator=1e6, decimal_places=4, verbose=True
)
df
df.to_csv('esa_area.csv')
Map = geemap.Map(center=[40, -100], zoom=4)
Map.add_basemap('HYBRID')

nlcd = ee.Image('USGS/NLCD_RELEASES/2019_REL/NLCD/2019')
landcover = nlcd.select('landcover')

Map.addLayer(landcover, {}, 'NLCD Land Cover 2019')
Map.add_legend(
    title="NLCD Land Cover Classification", builtin_legend='NLCD', height='465px'
)
Map
df = geemap.image_area_by_group(
    landcover, scale=1000, denominator=1e6, decimal_places=4, verbose=True
)
df
df.to_csv('nlcd_area.csv')

6.18. 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.19. WhiteboxTools#

pip install geemap[lidar]

6.19.1. Import libraries#

import os
import geemap
import whitebox

6.19.2. Set up whitebox#

wbt = whitebox.WhiteboxTools()
wbt.set_working_dir(os.getcwd())
wbt.set_verbose_mode(False)

6.19.3. Download sample data#

url = 'https://github.com/giswqs/data/raw/main/lidar/madison.laz'
if not os.path.exists('madison.laz'):
    geemap.download_file(url)

6.19.4. Read LAS/LAZ data#

laz = geemap.read_lidar('madison.laz')
laz
str(laz.header.version)

6.19.5. Upgrade file version#

las = geemap.convert_lidar(laz, file_version='1.4')
str(las.header.version)

6.19.6. Write LAS data#

geemap.write_lidar(las, 'madison.las')

6.19.7. Histogram analysis#

wbt.lidar_histogram('madison.las', 'histogram.html')

6.19.8. Visualize LiDAR data#

geemap.view_lidar('madison.las')

6.19.9. Remove outliers#

wbt.lidar_elevation_slice("madison.las", "madison_rm.las", minz=0, maxz=450)

6.19.10. Visualize LiDAR data after removing outliers#

geemap.view_lidar('madison_rm.las', cmap='terrain')

6.19.11. Create DSM#

wbt.lidar_digital_surface_model(
    'madison_rm.las', 'dsm.tif', resolution=1.0, minz=0, maxz=450
)
geemap.add_crs("dsm.tif", epsg=2255)

6.19.12. Visualize DSM#

Map = geemap.Map()
Map.add_raster('dsm.tif', palette='terrain', layer_name='DSM')
Map

6.19.13. Create DEM#

wbt.remove_off_terrain_objects('dsm.tif', 'dem.tif', filter=25, slope=15.0)

6.19.14. Visualize DEM#

Map.add_raster('dem.tif', palette='terrain', layer_name='DEM')
Map

6.19.15. Create CHM#

chm = wbt.subtract('dsm.tif', 'dem.tif', 'chm.tif')
Map.add_raster('chm.tif', palette='gist_earth', layer_name='CHM')
Map

6.20. Summary#

6.21. References#