Crop Mapping 2022

Cropland mapping with Google Earth Engine and geemap

Useful links

• Google Earth Engine
• Sign up for an Earth Engine account
• Awesome GEE Community Datasets
• Geemap website
• Geemap book
• Geemap YouTube videos

In [ ]:

# !pip install geemap

# !pip install geemap

Import libraries¶

In [ ]:

import ee
import geemap

import ee import geemap

ESA WordCover¶

The European Space Agency (ESA) WorldCover 10 m 2020 product provides a global land cover map for 2020 at 10 m resolution based on Sentinel-1 and Sentinel-2 data. The WorldCover product comes with 11 land cover classes and has been generated in the framework of the ESA WorldCover project, part of the 5th Earth Observation Envelope Programme (EOEP-5) of the European Space Agency.

• ESA WorldCover website
• EAS WroldCover in the Earth Engine Data Catalog
• User Manual and Validation Report

Using Web Map Services¶

The ESA WorldCover product can also be used within other websites or GIS clients by 'Web Map Services'. These services provide a direct link to the cached images and are the best option if you simply want to map the data and produce cartographic products. They are not suitable for analysis as the data are represented only as RGB images.

• WMTS: https://services.terrascope.be/wmts/v2
• WMS: https://services.terrascope.be/wms/v2
• Layers: WORLDCOVER_2020_MAP, WORLDCOVER_2020_S2_FCC, WORLDCOVER_2020_S2_TCC

In [ ]:

Map = geemap.Map()

esa_wms = "https://services.terrascope.be/wms/v2"  # The WMS URL
tcc_layer = "WORLDCOVER_2020_S2_TCC"  # The true color composite imagery
fcc_layer = "WORLDCOVER_2020_S2_FCC"  # The false color composite imagery
map_layer = "WORLDCOVER_2020_MAP"  # The land cover classification map

Map.add_wms_layer(esa_wms, layers=tcc_layer, name="True Color", attribution="ESA")
Map.add_wms_layer(esa_wms, layers=fcc_layer, name="False Color", attribution="ESA")
Map.add_wms_layer(esa_wms, layers=map_layer, name="Classification", attribution="ESA")

Map.add_legend(title="ESA Land Cover", builtin_legend="ESA_WorldCover")
Map

Map = geemap.Map() esa_wms = "https://services.terrascope.be/wms/v2" # The WMS URL tcc_layer = "WORLDCOVER_2020_S2_TCC" # The true color composite imagery fcc_layer = "WORLDCOVER_2020_S2_FCC" # The false color composite imagery map_layer = "WORLDCOVER_2020_MAP" # The land cover classification map Map.add_wms_layer(esa_wms, layers=tcc_layer, name="True Color", attribution="ESA") Map.add_wms_layer(esa_wms, layers=fcc_layer, name="False Color", attribution="ESA") Map.add_wms_layer(esa_wms, layers=map_layer, name="Classification", attribution="ESA") Map.add_legend(title="ESA Land Cover", builtin_legend="ESA_WorldCover") Map

Using Earth Engine¶

• EAS WroldCover in the Earth Engine Data Catalog

In [ ]:

Map = geemap.Map()
Map.add_basemap("HYBRID")

esa = ee.ImageCollection("ESA/WorldCover/v100").first()
esa_vis = {"bands": ["Map"]}

Map.addLayer(esa, esa_vis, "ESA Land Cover")
Map.add_legend(title="ESA Land Cover", builtin_legend="ESA_WorldCover")

Map

Map = geemap.Map() Map.add_basemap("HYBRID") esa = ee.ImageCollection("ESA/WorldCover/v100").first() esa_vis = {"bands": ["Map"]} Map.addLayer(esa, esa_vis, "ESA Land Cover") Map.add_legend(title="ESA Land Cover", builtin_legend="ESA_WorldCover") Map

Creating charts¶

In [ ]:

histogram = geemap.image_histogram(
    esa, scale=1000, x_label="Land Cover Type", y_label="Area (km2)"
)
histogram

histogram = geemap.image_histogram( esa, scale=1000, x_label="Land Cover Type", y_label="Area (km2)" ) histogram

In [ ]:

df = geemap.image_histogram(esa, scale=1000, return_df=True)
df

df = geemap.image_histogram(esa, scale=1000, return_df=True) df

In [ ]:

esa_labels = list(geemap.builtin_legends["ESA_WorldCover"].keys())
esa_labels

esa_labels = list(geemap.builtin_legends["ESA_WorldCover"].keys()) esa_labels

In [ ]:

df["label"] = esa_labels
df

df["label"] = esa_labels df

In [ ]:

round(df["value"].sum() / 1e6, 2)

round(df["value"].sum() / 1e6, 2)

In [ ]:

geemap.bar_chart(
    df, x="label", y="value", x_label="Land Cover Type", y_label="Area (km2)"
)

geemap.bar_chart( df, x="label", y="value", x_label="Land Cover Type", y_label="Area (km2)" )

In [ ]:

geemap.pie_chart(df, names="label", values="value", height=500)

geemap.pie_chart(df, names="label", values="value", height=500)

Adding Administrative Boundaries¶

In [ ]:

countries = ee.FeatureCollection(geemap.examples.get_ee_path("countries"))
africa = countries.filter(ee.Filter.eq("CONTINENT", "Africa"))
style = {"fillColor": "00000000"}
Map.addLayer(countries.style(**style), {}, "Countries", False)
Map.addLayer(africa.style(**style), {}, "Africa")
Map.centerObject(africa)
Map

countries = ee.FeatureCollection(geemap.examples.get_ee_path("countries")) africa = countries.filter(ee.Filter.eq("CONTINENT", "Africa")) style = {"fillColor": "00000000"} Map.addLayer(countries.style(**style), {}, "Countries", False) Map.addLayer(africa.style(**style), {}, "Africa") Map.centerObject(africa) Map

Extracting Croplands¶

In [ ]:

cropland = esa.eq(40).clipToCollection(africa).selfMask()
Map.addLayer(cropland, {"palette": ["f096ff"]}, "Cropland")
Map.show_layer(name="ESA Land Cover", show=False)

cropland = esa.eq(40).clipToCollection(africa).selfMask() Map.addLayer(cropland, {"palette": ["f096ff"]}, "Cropland") Map.show_layer(name="ESA Land Cover", show=False)

Zonal Statistics¶

In [ ]:

geemap.zonal_stats(cropland, africa, "esa_cropland.csv", stat_type="SUM", scale=1000)

geemap.zonal_stats(cropland, africa, "esa_cropland.csv", stat_type="SUM", scale=1000)

In [ ]:

df = geemap.csv_to_df("esa_cropland.csv")
df.head()

df = geemap.csv_to_df("esa_cropland.csv") df.head()

In [ ]:

geemap.bar_chart(
    df, x="NAME", y="sum", max_rows=30, x_label="Country", y_label="Area (km2)"
)

geemap.bar_chart( df, x="NAME", y="sum", max_rows=30, x_label="Country", y_label="Area (km2)" )

In [ ]:

geemap.pie_chart(df, names="NAME", values="sum", max_rows=20, height=500)

geemap.pie_chart(df, names="NAME", values="sum", max_rows=20, height=500)

ESRI GLobal Land Cover¶

The ESRI GLobal Land Cover dataset is a global map of land use/land cover (LULC) derived from ESA Sentinel-2 imagery at 10m resolution. Each year is generated from Impact Observatory’s deep learning AI land classification model used a massive training dataset of billions of human-labeled image pixels developed by the National Geographic Society. The global maps were produced by applying this model to the Sentinel-2 scene collection on Microsoft’s Planetary Computer, processing over 400,000 Earth observations per year.

• https://livingatlas.arcgis.com/landcover/
• https://www.arcgis.com/home/item.html?id=d3da5dd386d140cf93fc9ecbf8da5e31
• https://samapriya.github.io/awesome-gee-community-datasets/projects/S2TSLULC/

Using Awesome GEE Community Datasets¶

In [ ]:

Map = geemap.Map()
Map.add_basemap("HYBRID")

esri = ee.ImageCollection(
    "projects/sat-io/open-datasets/landcover/ESRI_Global-LULC_10m_TS"
)

esri_2017 = esri.filterDate("2017-01-01", "2017-12-31").mosaic()
esri_2018 = esri.filterDate("2018-01-01", "2018-12-31").mosaic()
esri_2019 = esri.filterDate("2019-01-01", "2019-12-31").mosaic()
esri_2020 = esri.filterDate("2020-01-01", "2020-12-31").mosaic()
esri_2021 = esri.filterDate("2021-01-01", "2021-12-31").mosaic()

esri_vis = {"min": 1, "max": 11, "palette": "esri_lulc"}

Map.addLayer(esri_2017, esri_vis, "ESRI LULC 2017")
Map.addLayer(esri_2018, esri_vis, "ESRI LULC 2018")
Map.addLayer(esri_2019, esri_vis, "ESRI LULC 2019")
Map.addLayer(esri_2020, esri_vis, "ESRI LULC 2020")
Map.addLayer(esri_2021, esri_vis, "ESRI LULC 2021")

Map.add_legend(title="ESRI Land Cover", builtin_legend="ESRI_LandCover_TS")
Map

Map = geemap.Map() Map.add_basemap("HYBRID") esri = ee.ImageCollection( "projects/sat-io/open-datasets/landcover/ESRI_Global-LULC_10m_TS" ) esri_2017 = esri.filterDate("2017-01-01", "2017-12-31").mosaic() esri_2018 = esri.filterDate("2018-01-01", "2018-12-31").mosaic() esri_2019 = esri.filterDate("2019-01-01", "2019-12-31").mosaic() esri_2020 = esri.filterDate("2020-01-01", "2020-12-31").mosaic() esri_2021 = esri.filterDate("2021-01-01", "2021-12-31").mosaic() esri_vis = {"min": 1, "max": 11, "palette": "esri_lulc"} Map.addLayer(esri_2017, esri_vis, "ESRI LULC 2017") Map.addLayer(esri_2018, esri_vis, "ESRI LULC 2018") Map.addLayer(esri_2019, esri_vis, "ESRI LULC 2019") Map.addLayer(esri_2020, esri_vis, "ESRI LULC 2020") Map.addLayer(esri_2021, esri_vis, "ESRI LULC 2021") Map.add_legend(title="ESRI Land Cover", builtin_legend="ESRI_LandCover_TS") Map

Using Timeseries Inspector¶

In [ ]:

images = ee.List([esri_2017, esri_2018, esri_2019, esri_2020, esri_2021])
collection = ee.ImageCollection.fromImages(images)

images = ee.List([esri_2017, esri_2018, esri_2019, esri_2020, esri_2021]) collection = ee.ImageCollection.fromImages(images)

In [ ]:

years = [str(year) for year in range(2017, 2022)]
years

years = [str(year) for year in range(2017, 2022)] years

In [ ]:

Map = geemap.Map()
Map.ts_inspector(collection, years, esri_vis, width="80px")
Map.add_legend(title="ESRI Land Cover", builtin_legend="ESRI_LandCover_TS")
Map

Map = geemap.Map() Map.ts_inspector(collection, years, esri_vis, width="80px") Map.add_legend(title="ESRI Land Cover", builtin_legend="ESRI_LandCover_TS") Map

Extracting Croplands¶

In [ ]:

countries = ee.FeatureCollection(geemap.examples.get_ee_path("countries"))
africa = countries.filter(ee.Filter.eq("CONTINENT", "Africa"))

countries = ee.FeatureCollection(geemap.examples.get_ee_path("countries")) africa = countries.filter(ee.Filter.eq("CONTINENT", "Africa"))

In [ ]:

cropland_col = collection.map(lambda img: img.eq(5).clipToCollection(africa).selfMask())
cropland_ts = cropland_col.toBands().rename(years)

cropland_col = collection.map(lambda img: img.eq(5).clipToCollection(africa).selfMask()) cropland_ts = cropland_col.toBands().rename(years)

In [ ]:

Map = geemap.Map()

style = {"fillColor": "00000000"}
Map.addLayer(cropland_col.first(), {"palette": ["#ab6c28"]}, "first")
Map.addLayer(countries.style(**style), {}, "Countries", False)
Map.addLayer(africa.style(**style), {}, "Africa")
Map.centerObject(africa)

Map

Map = geemap.Map() style = {"fillColor": "00000000"} Map.addLayer(cropland_col.first(), {"palette": ["#ab6c28"]}, "first") Map.addLayer(countries.style(**style), {}, "Countries", False) Map.addLayer(africa.style(**style), {}, "Africa") Map.centerObject(africa) Map

In [ ]:

cropland_ts.bandNames().getInfo()

cropland_ts.bandNames().getInfo()

Zonal Statistics¶

In [ ]:

geemap.zonal_stats(
    cropland_ts, africa, "esri_cropland.csv", stat_type="SUM", scale=1000
)

geemap.zonal_stats( cropland_ts, africa, "esri_cropland.csv", stat_type="SUM", scale=1000 )

In [ ]:

df = geemap.csv_to_df("esri_cropland.csv")
df.head()

df = geemap.csv_to_df("esri_cropland.csv") df.head()

In [ ]:

geemap.bar_chart(df, x="NAME", y=years, max_rows=20, legend_title="Year")

geemap.bar_chart(df, x="NAME", y=years, max_rows=20, legend_title="Year")

In [ ]:

geemap.pie_chart(df, names="NAME", values="2020", max_rows=20, height=500)

geemap.pie_chart(df, names="NAME", values="2020", max_rows=20, height=500)

Analyzing Cropland Gain and Loss¶

In [ ]:

Map = geemap.Map()
Map.add_basemap("HYBRID")

cropland_2017 = esri_2017.eq(5).selfMask()
cropland_2021 = esri_2021.eq(5).selfMask()

cropland_gain = esri_2017.neq(5).And(esri_2021.eq(5)).selfMask()
cropland_loss = esri_2017.eq(5).And(esri_2021.neq(5)).selfMask()

Map.addLayer(cropland_2017, {"palette": "brown"}, "Cropland 2017", False)
Map.addLayer(cropland_2021, {"palette": "cyan"}, "Cropland 2021", False)

Map.addLayer(cropland_gain, {"palette": "yellow"}, "Cropland gain")
Map.addLayer(cropland_loss, {"palette": "red"}, "Cropland loss")
Map

Map = geemap.Map() Map.add_basemap("HYBRID") cropland_2017 = esri_2017.eq(5).selfMask() cropland_2021 = esri_2021.eq(5).selfMask() cropland_gain = esri_2017.neq(5).And(esri_2021.eq(5)).selfMask() cropland_loss = esri_2017.eq(5).And(esri_2021.neq(5)).selfMask() Map.addLayer(cropland_2017, {"palette": "brown"}, "Cropland 2017", False) Map.addLayer(cropland_2021, {"palette": "cyan"}, "Cropland 2021", False) Map.addLayer(cropland_gain, {"palette": "yellow"}, "Cropland gain") Map.addLayer(cropland_loss, {"palette": "red"}, "Cropland loss") Map

In [ ]:

geemap.zonal_stats(
    cropland_gain,
    countries,
    "esri_cropland_gain.csv",
    stat_type="SUM",
    scale=1000,
)

geemap.zonal_stats( cropland_gain, countries, "esri_cropland_gain.csv", stat_type="SUM", scale=1000, )

In [ ]:

df = geemap.csv_to_df("esri_cropland_gain.csv")
df.head()

df = geemap.csv_to_df("esri_cropland_gain.csv") df.head()

In [ ]:

geemap.bar_chart(
    df,
    x="NAME",
    y="sum",
    max_rows=30,
    x_label="Country",
    y_label="Area (km2)",
    title="Cropland Gain",
)

geemap.bar_chart( df, x="NAME", y="sum", max_rows=30, x_label="Country", y_label="Area (km2)", title="Cropland Gain", )

In [ ]:

geemap.pie_chart(
    df, names="NAME", values="sum", max_rows=30, height=500, title="Cropland Gain"
)

geemap.pie_chart( df, names="NAME", values="sum", max_rows=30, height=500, title="Cropland Gain" )

In [ ]:

geemap.zonal_stats(
    cropland_loss,
    countries,
    "esri_cropland_loss.csv",
    stat_type="SUM",
    scale=1000,
)

geemap.zonal_stats( cropland_loss, countries, "esri_cropland_loss.csv", stat_type="SUM", scale=1000, )

In [ ]:

df = geemap.csv_to_df("esri_cropland_loss.csv")
df.head()

df = geemap.csv_to_df("esri_cropland_loss.csv") df.head()

In [ ]:

geemap.bar_chart(
    df,
    x="NAME",
    y="sum",
    max_rows=30,
    x_label="Country",
    y_label="Area (km2)",
    title="Cropland Loss",
)

geemap.bar_chart( df, x="NAME", y="sum", max_rows=30, x_label="Country", y_label="Area (km2)", title="Cropland Loss", )

In [ ]:

geemap.pie_chart(
    df, names="NAME", values="sum", max_rows=30, height=500, title="Cropland Loss"
)

geemap.pie_chart( df, names="NAME", values="sum", max_rows=30, height=500, title="Cropland Loss" )

Dynamic World Land Cover¶

Dynamic World is a near realtime 10m resolution global land use land cover dataset, produced using deep learning, freely available and openly licensed. As a result of leveraging a novel deep learning approach, based on Sentinel-2 Top of Atmosphere, Dynamic World offers global land cover updating every 2-5 days depending on location.

• Dynamic World Website
• Dynamic World datasets on Earth Engine

Classification and Probability¶

In [ ]:

Map = geemap.Map()

region = ee.Geometry.BBox(-179, -89, 179, 89)
start_date = "2021-01-01"
end_date = "2022-01-01"

dw_class = geemap.dynamic_world(region, start_date, end_date, return_type="class")
dw = geemap.dynamic_world(region, start_date, end_date, return_type="hillshade")

dw_vis = {"min": 0, "max": 8, "palette": "dw"}

Map.addLayer(dw_class, dw_vis, "DW Land Cover", False)
Map.addLayer(dw, {}, "DW Land Cover Hillshade")

Map.add_legend(title="Dynamic World Land Cover", builtin_legend="Dynamic_World")
Map.setCenter(-88.9088, 43.0006, 12)
Map

Map = geemap.Map() region = ee.Geometry.BBox(-179, -89, 179, 89) start_date = "2021-01-01" end_date = "2022-01-01" dw_class = geemap.dynamic_world(region, start_date, end_date, return_type="class") dw = geemap.dynamic_world(region, start_date, end_date, return_type="hillshade") dw_vis = {"min": 0, "max": 8, "palette": "dw"} Map.addLayer(dw_class, dw_vis, "DW Land Cover", False) Map.addLayer(dw, {}, "DW Land Cover Hillshade") Map.add_legend(title="Dynamic World Land Cover", builtin_legend="Dynamic_World") Map.setCenter(-88.9088, 43.0006, 12) Map

ESA Land Cover vs. Dynamic World¶

In [ ]:

Map = geemap.Map(center=[39.3322, -106.7349], zoom=10)

left_layer = geemap.ee_tile_layer(esa, esa_vis, "ESA Land Cover")
right_layer = geemap.ee_tile_layer(dw, {}, "Dynamic World Land Cover")

Map.split_map(left_layer, right_layer)
Map.add_legend(
    title="ESA Land Cover", builtin_legend="ESA_WorldCover", position="bottomleft"
)
Map.add_legend(
    title="Dynamic World Land Cover",
    builtin_legend="Dynamic_World",
    position="bottomright",
)
Map.setCenter(-88.9088, 43.0006, 12)

Map

Map = geemap.Map(center=[39.3322, -106.7349], zoom=10) left_layer = geemap.ee_tile_layer(esa, esa_vis, "ESA Land Cover") right_layer = geemap.ee_tile_layer(dw, {}, "Dynamic World Land Cover") Map.split_map(left_layer, right_layer) Map.add_legend( title="ESA Land Cover", builtin_legend="ESA_WorldCover", position="bottomleft" ) Map.add_legend( title="Dynamic World Land Cover", builtin_legend="Dynamic_World", position="bottomright", ) Map.setCenter(-88.9088, 43.0006, 12) Map

ESRI Land Cover vs. Dynamic World¶

In [ ]:

Map = geemap.Map(center=[-89.3998, 43.0886], zoom=10)

left_layer = geemap.ee_tile_layer(esri_2021, esri_vis, "ESRI Land Cover")
right_layer = geemap.ee_tile_layer(dw, {}, "Dynamic World Land Cover")

Map.split_map(left_layer, right_layer)
Map.add_legend(
    title="ESRI Land Cover", builtin_legend="ESRI_LandCover", position="bottomleft"
)
Map.add_legend(
    title="Dynamic World Land Cover",
    builtin_legend="Dynamic_World",
    position="bottomright",
)
Map.setCenter(-88.9088, 43.0006, 12)

Map

Map = geemap.Map(center=[-89.3998, 43.0886], zoom=10) left_layer = geemap.ee_tile_layer(esri_2021, esri_vis, "ESRI Land Cover") right_layer = geemap.ee_tile_layer(dw, {}, "Dynamic World Land Cover") Map.split_map(left_layer, right_layer) Map.add_legend( title="ESRI Land Cover", builtin_legend="ESRI_LandCover", position="bottomleft" ) Map.add_legend( title="Dynamic World Land Cover", builtin_legend="Dynamic_World", position="bottomright", ) Map.setCenter(-88.9088, 43.0006, 12) Map