Change in Earth’s albedo and its dependence on land cover changes in the past 20 years¶

Content creators: Oz Kira, Julius Bamah

Content reviewers: Yuhan Douglas Rao, Abigail Bodner

Content editors: Zane Mitrevica, Natalie Steinemann, Jenna Pearson, Chi Zhang, Ohad Zivan

Production editors: Wesley Banfield, Jenna Pearson, Chi Zhang, Ohad Zivan

Our 2023 Sponsors: NASA TOPS and Google DeepMind

# @title #**Project background** 
#This will be a short video introducing the content creator(s) and motivating the research direction of the template.
#The Tech team will add code to format and display the video

The global radiative budget is affected by land cover. Regarding vegetation land cover (e.g., forests, grasslands, agricultural fields, etc.), vegetation sequesters carbon, which reduces the greenhouse effect but absorbs more radiation and reduces earth albedo, which counteracts carbon sequestration.

In this project, we will evaluate the albedo change vs. carbon sequestration over the past years. In addition, we will track significant land use changes, specifically the creation and abandonment of agricultural land.

In this project, you will have the opportunity to explore terrestrial remote sensing (recall our W1D3 tutorial on remote sensing) and meteorological data from GLASS and ERA5. The datasets will provide information on reflectance, albedo, meteorological variables, and land cover changes in your region of interest. We encourage you to investigate the relationships between these variables and their impact on the global radiative budget. Moreover, you can track agricultural land abandonment and analyze its potential connection to climate change. This project aligns well with the topics covered in W2D3, which you are encouraged to explore further.

Project template¶

template

Data exploration notebook¶

Project setup¶

Please run the following cells to install the necessary libarries into your Jupyter notebook!

# google colab installs
# !pip install cartopy
# !pip install DateTime 
# !pip install matplotlib
# !pip install pyhdf
# !pip install numpy
# !pip install pandas
# !pip install modis-tools

# the further information on the MODIS data can be found here : 
#Import the libraries
import numpy as np
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pooch
import xarray as xr
import os

Global Land Surface Satellite (GLASS) Dataset¶

The Global Land Surface Satellite (GLASS) datasets primarily based on NASA’s Advanced Very High Resolution Radiometer (AVHRR) long-term data record (LTDR) and Moderate Resolution Imaging Spectroradiometer (MODIS) data, in conjunction with other satellite data and ancillary information.

Currently, there are more than dozens of GLASS products are officially released, including leaf area index, fraction of green vegetation coverage, gross primary production, broadband albedo, land surface temperature, evapotranspiration, and so on.

Here we provide you the datasets of GLASS from 1982 to 2015, a 34-year long annual dynamics of global land cover (GLASS-GLC) at 5 km resolution. In this datasets, there are 7 classes, including cropland, forest, grassland, shrubland, tundra, barren land, and snow/ice. The annual global land cover map (5 km) is presented in a GeoTIFF file format named in the form of ‘GLASS-GLC_7classes_year’ with a WGS 84 projection. The relationship between the labels in the files and the 7 land cover classes is shown in the following table

You can refer to this paper for detailed description of this.ts

# Table 1 Classification system, with 7 land cover classes. From paper https://www.earth-syst-sci-data-discuss.net/essd-2019-23
import pandas as pd
from IPython.display import display, HTML, Markdown
# Data as list of dictionaries
classification_system = [
    {"Label": 10, "Class": "Cropland", "Subclass": "Rice paddy", "Description": ""},
    {"Label": 10, "Class": "Cropland", "Subclass": "Greenhouse", "Description": ""},
    {"Label": 10, "Class": "Cropland", "Subclass": "Other farmland", "Description": ""},
    {"Label": 10, "Class": "Cropland", "Subclass": "Orchard", "Description": ""},
    {"Label": 10, "Class": "Cropland", "Subclass": "Bare farmland", "Description": ""},
    {"Label": 20, "Class": "Forest", "Subclass": "Broadleaf, leaf-on", "Description": "Tree cover≥10%; Height>5m; For mixed leaf, neither coniferous nor broadleaf types exceed 60%"},
    {"Label": 20, "Class": "Forest", "Subclass": "Broadleaf, leaf-off", "Description": ""},
    {"Label": 20, "Class": "Forest", "Subclass": "Needle-leaf, leaf-on", "Description": ""},
    {"Label": 20, "Class": "Forest", "Subclass": "Needle-leaf, leaf-off", "Description": ""},
    {"Label": 20, "Class": "Forest", "Subclass": "Mixed leaf type, leaf-on", "Description": ""},
    {"Label": 20, "Class": "Forest", "Subclass": "Mixed leaf type, leaf-off", "Description": ""},
    {"Label": 30, "Class": "Grassland", "Subclass": "Pasture, leaf-on", "Description": "Canopy cover≥20%"},
    {"Label": 30, "Class": "Grassland", "Subclass": "Natural grassland, leaf-on", "Description": ""},
    {"Label": 30, "Class": "Grassland", "Subclass": "Grassland, leaf-off", "Description": ""},
    {"Label": 40, "Class": "Shrubland", "Subclass": "Shrub cover, leaf-on", "Description": "Canopy cover≥20%; Height<5m"},
    {"Label": 40, "Class": "Shrubland", "Subclass": "Shrub cover, leaf-off", "Description": ""},
    {"Label": 70, "Class": "Tundra", "Subclass": "Shrub and brush tundra", "Description": ""},
    {"Label": 70, "Class": "Tundra", "Subclass": "Herbaceous tundra", "Description": ""},
    {"Label": 90, "Class": "Barren land", "Subclass": "Barren land", "Description": "Vegetation cover<10%"},
    {"Label": 100, "Class": "Snow/Ice", "Subclass": "Snow", "Description": ""},
    {"Label": 100, "Class": "Snow/Ice", "Subclass": "Ice", "Description": ""},
    {"Label": 0, "Class": "No data", "Subclass": "", "Description": ""}
]

df = pd.DataFrame(classification_system)
pd.set_option('display.max_colwidth', None)
html = df.to_html(index=False)
title_md = "### Table 1 GLASS classification system with 7 land cover classes. From [this paper](https://www.earth-syst-sci-data-discuss.net/essd-2019-23)."
display(Markdown(title_md))
display(HTML(html))

System Message: WARNING/2 (/home/yuda/code/climatematch/course-content/book/projects/project-notebooks/Surface_albedo_and_land_use.ipynb, line 1)

Non-consecutive header level increase; 0 to 3 [myst.header]

Label	Class	Subclass	Description
10	Cropland	Rice paddy
10	Cropland	Greenhouse
10	Cropland	Other farmland
10	Cropland	Orchard
10	Cropland	Bare farmland
20	Forest	Broadleaf, leaf-on	Tree cover≥10%; Height>5m; For mixed leaf, neither coniferous nor broadleaf types exceed 60%
20	Forest	Broadleaf, leaf-off
20	Forest	Needle-leaf, leaf-on
20	Forest	Needle-leaf, leaf-off
20	Forest	Mixed leaf type, leaf-on
20	Forest	Mixed leaf type, leaf-off
30	Grassland	Pasture, leaf-on	Canopy cover≥20%
30	Grassland	Natural grassland, leaf-on
30	Grassland	Grassland, leaf-off
40	Shrubland	Shrub cover, leaf-on	Canopy cover≥20%; Height<5m
40	Shrubland	Shrub cover, leaf-off
70	Tundra	Shrub and brush tundra
70	Tundra	Herbaceous tundra
90	Barren land	Barren land	Vegetation cover<10%
100	Snow/Ice	Snow
100	Snow/Ice	Ice
0	No data

Alternative Approach to Obtain Land Cover Data and Net Primary Production (NPP) Data from MODIS¶

MODIS (Moderate Resolution Imaging Spectroradiometer) is a key instrument aboard the Terra (originally known as EOS AM-1) and Aqua (originally known as EOS PM-1) satellites. Terra’s orbit around the Earth is timed so that it passes from north to south across the equator in the morning, while Aqua passes south to north over the equator in the afternoon. Terra MODIS and Aqua MODIS are viewing the entire Earth’s surface every 1 to 2 days, acquiring data in 36 spectral bands, or groups of wavelengths (see MODIS Technical Specifications). These data will improve our understanding of global dynamics and processes occurring on the land, in the oceans, and in the lower atmosphere. MODIS is playing a vital role in the development of validated, global, interactive Earth system models able to predict global change accurately enough to assist policy makers in making sound decisions concerning the protection of our environment (https://modis.gsfc.nasa.gov/data/).

The following procedures must be ensured:

Register an account on register at NASA earth data portal: https://urs.earthdata.nasa.gov/users/new
pip install the following libraries into your Jupyter notebook

The following links could help you in preprocessing the MODIS dataset

https://www.earthdatascience.org/courses/use-data-open-source-python/multispectral-remote-sensing/modis-data-in-python/
http://www.hdfeos.org/zoo/LAADS/MYD08_D3.A2009001.006.2015088160849.hdf.py
https://www.moonbooks.org/Articles/How-to-read-a-MODIS-HDF-file-using-python-/

Before running the code please read this.

The following papers could assist in processing of MODIS dataset:

https://www.mdpi.com/2072-4292/8/7/554
https://www.tandfonline.com/doi/full/10.1080/01431161.2018.1430913
https://www.mdpi.com/2072-4292/6/6/5368

## Alternative Approach to Obtain Land Cover Data and Net Primary Production (NPP) Data from MODIS
# #Import the libraries

# from modis_tools.auth import ModisSession
# from modis_tools.resources import CollectionApi, GranuleApi
# from modis_tools.granule_handler import GranuleHandler
# from typing_extensions import Literal

## Alternative Approach to Obtain Land Cover Data and Net Primary Production (NPP) Data from MODIS
# #download_modis.py
# username = ""  # Update this line with your username 
# password = ""  # Update this line with your password
# # 1) connect to earthdata
# session = ModisSession(username=username, password=password)

# # 2) Query the MODIS catalog for collections
# collection_client = CollectionApi(session=session)
# collections = collection_client.query(short_name="", version="") # Collection short name + version
# # Query the selected collection for granules
# granule_client = GranuleApi.from_collection(collections[0], session=session)

# # 3) Filter the selected granules via spatial and temporal parameters
# Israel_bbox = [] # format [x_min, y_min, x_max, y_max] # add min lon,min lat, max lon and max_lat input your preferred location
# Israel_granules = granule_client.query(start_date="", end_date="", bounding_box=Israel_bbox)   #choose the start and end dates for the year-month-day you prefer

# # 4) Download the granules
# GranuleHandler.download_from_granules(Israel_granules, session, threads=-1) 
# #NB the file format will be downloaded in hdf 

Gross Primary Production (GPP) Datasets¶

In the tutorial, you learned about Net primary production (NPP). Another similar key ecosystem process is Gross Primary Production (GPP), which is the total amount of carbon compounds produced by photosynthesis of plants in an ecosystem in a given period of time. NPP is equal to GPP minus energy used by primary producers for respiration. GPP is the amount of energy from light converted to chemical energy per unit of time, while NPP is a measure of the net CO2 sequestered by vegetation.

Here we provide you the GPP datasets VODCA2GPP via here and this paper. This dataset estimates GPP based on VOD (Microwave Vegetation Optical Depth (VOD), a measure of the attenuation of microwave radiation caused by vegetation and thus relates to the total vegetation water content) using Generalized Additive Models. This VODCA2GPP is based on the period from 2007 to 2015 and uses VOD data from C-, X- and Ku-band and various GPP data sets. The data sets have different temporal coverage, which is summarized for VOD and GPP data in Table 2.

import pandas as pd
from IPython.display import display, HTML, Markdown

# Data as list of dictionaries
data = [
    {"Sensor": "AMSR-E", "Time period used": "Jun 2002–Oct 2011", "AECT": "13:30", "C-band [GHz]": "6.93", "X-band [GHz]": "10.65", "Ku-band [GHz]": "18.7", "Reference": "Van der Schalie et al. (2017)"},
    {"Sensor": "AMSR2", "Time period used": "Jul 2012–Dec 2020", "AECT": "13:30", "C-band [GHz]": "6.93, 7.30", "X-band [GHz]": "10.65", "Ku-band [GHz]": "18.7", "Reference": "Van der Schalie et al. (2017)"},
    {"Sensor": "AMSR2", "Time period used": "Jul 2012–Aug 2017 (Ku-band)", "AECT": "", "C-band [GHz]": "", "X-band [GHz]": "", "Ku-band [GHz]": "", "Reference": ""},
    {"Sensor": "SSM/I F08", "Time period used": "Jul 1987–Dec 1991", "AECT": "18:15", "C-band [GHz]": "", "X-band [GHz]": "", "Ku-band [GHz]": "19.35", "Reference": "Owe et al. (2008)"},
    {"Sensor": "SSM/I F11", "Time period used": "Dec 1991–May 1995", "AECT": "17:00–18:15", "C-band [GHz]": "", "X-band [GHz]": "", "Ku-band [GHz]": "19.35", "Reference": "Owe et al. (2008)"},
    {"Sensor": "SSM/I F13", "Time period used": "May 1995–Apr 2009", "AECT": "17:45–18:40", "C-band [GHz]": "", "X-band [GHz]": "", "Ku-band [GHz]": "19.35", "Reference": "Owe et al. (2008)"},
    {"Sensor": "TMI", "Time period used": "Dec 1997–Apr 2015", "AECT": "Asynchronous", "C-band [GHz]": "", "X-band [GHz]": "10.65", "Ku-band [GHz]": "19.35", "Reference": "Owe et al. (2008), Van der Schalie et al. (2017)"},
    {"Sensor": "WindSat", "Time period used": "Feb 2003–Jul 2012", "AECT": "18:00", "C-band [GHz]": "6.8", "X-band [GHz]": "10.7", "Ku-band [GHz]": "18.7", "Reference": "Owe et al. (2008), Van der Schalie et al. (2017)"}
]

# Convert list of dictionaries into a DataFrame
df = pd.DataFrame(data)

# Set max_colwidth to None so that the full content of Description column is displayed
pd.set_option('display.max_colwidth', None)

# Convert the DataFrame to HTML and hide the index
html = df.to_html(index=False)

# Table title and link to paper in Markdown
title_md = "### Table 2: Variable names, data sets, and additional information of VODCA2GPP datasets. [Reference](https://essd.copernicus.org/articles/14/1063/2022/#section6)"

# Display the title
display(Markdown(title_md))

# Display the table
display(HTML(html))

System Message: WARNING/2 (/home/yuda/code/climatematch/course-content/book/projects/project-notebooks/Surface_albedo_and_land_use.ipynb, line 1)

Non-consecutive header level increase; 0 to 3 [myst.header]

Sensor	Time period used	AECT	C-band [GHz]	X-band [GHz]	Ku-band [GHz]	Reference
AMSR-E	Jun 2002–Oct 2011	13:30	6.93	10.65	18.7	Van der Schalie et al. (2017)
AMSR2	Jul 2012–Dec 2020	13:30	6.93, 7.30	10.65	18.7	Van der Schalie et al. (2017)
AMSR2	Jul 2012–Aug 2017 (Ku-band)
SSM/I F08	Jul 1987–Dec 1991	18:15			19.35	Owe et al. (2008)
SSM/I F11	Dec 1991–May 1995	17:00–18:15			19.35	Owe et al. (2008)
SSM/I F13	May 1995–Apr 2009	17:45–18:40			19.35	Owe et al. (2008)
TMI	Dec 1997–Apr 2015	Asynchronous		10.65	19.35	Owe et al. (2008), Van der Schalie et al. (2017)
WindSat	Feb 2003–Jul 2012	18:00	6.8	10.7	18.7	Owe et al. (2008), Van der Schalie et al. (2017)

Net Primary Production (NPP) Datasets¶

Alternatively, you can access the NPP data by downloading the MCD12Q1.061 MODIS Land Cover Type Yearly Global 500m using the instruction provided earlier. The description of this dataset can be viewed here.

# gpp data is acquired from this work: https://essd.copernicus.org/articles/14/1063/2022/
url_GPP='https://researchdata.tuwien.ac.at/records/1k7aj-bdz35/files/VODCA2GPP_v1.nc?download=1'
ds_GPP=xr.open_dataset(pooch.retrieve(url_GPP,known_hash=None))
ds_GPP

Downloading data from 'https://researchdata.tuwien.ac.at/records/1k7aj-bdz35/files/VODCA2GPP_v1.nc?download=1' to file '/home/yuda/.cache/pooch/6e247026b602773579622b33be872c40-VODCA2GPP_v1.nc'.

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[8], line 3
      1 # gpp data is acquired from this work: https://essd.copernicus.org/articles/14/1063/2022/
      2 url_GPP='https://researchdata.tuwien.ac.at/records/1k7aj-bdz35/files/VODCA2GPP_v1.nc?download=1'
----> 3 ds_GPP=xr.open_dataset(pooch.retrieve(url_GPP,known_hash=None))
      4 ds_GPP

File ~/code/climatematch/course-content/venv/lib/python3.8/site-packages/pooch/core.py:239, in retrieve(url, known_hash, fname, path, processor, downloader, progressbar)
    236 if downloader is None:
    237     downloader = choose_downloader(url, progressbar=progressbar)
--> 239 stream_download(url, full_path, known_hash, downloader, pooch=None)
    241 if known_hash is None:
    242     get_logger().info(
    243         "SHA256 hash of downloaded file: %s\n"
    244         "Use this value as the 'known_hash' argument of 'pooch.retrieve'"
   (...)
    247         file_hash(str(full_path)),
    248     )

File ~/code/climatematch/course-content/venv/lib/python3.8/site-packages/pooch/core.py:803, in stream_download(url, fname, known_hash, downloader, pooch, retry_if_failed)
    799 try:
    800     # Stream the file to a temporary so that we can safely check its
    801     # hash before overwriting the original.
    802     with temporary_file(path=str(fname.parent)) as tmp:
--> 803         downloader(url, tmp, pooch)
    804         hash_matches(tmp, known_hash, strict=True, source=str(fname.name))
    805         shutil.move(tmp, str(fname))

File ~/code/climatematch/course-content/venv/lib/python3.8/site-packages/pooch/downloaders.py:226, in HTTPDownloader.__call__(self, url, output_file, pooch, check_only)
    224     progress = self.progressbar
    225     progress.total = total
--> 226 for chunk in content:
    227     if chunk:
    228         output_file.write(chunk)

File ~/code/climatematch/course-content/venv/lib/python3.8/site-packages/requests/models.py:816, in Response.iter_content.<locals>.generate()
    814 if hasattr(self.raw, "stream"):
    815     try:
--> 816         yield from self.raw.stream(chunk_size, decode_content=True)
    817     except ProtocolError as e:
    818         raise ChunkedEncodingError(e)

File ~/code/climatematch/course-content/venv/lib/python3.8/site-packages/urllib3/response.py:628, in HTTPResponse.stream(self, amt, decode_content)
    626 else:
    627     while not is_fp_closed(self._fp):
--> 628         data = self.read(amt=amt, decode_content=decode_content)
    630         if data:
    631             yield data

File ~/code/climatematch/course-content/venv/lib/python3.8/site-packages/urllib3/response.py:567, in HTTPResponse.read(self, amt, decode_content, cache_content)
    564 fp_closed = getattr(self._fp, "closed", False)
    566 with self._error_catcher():
--> 567     data = self._fp_read(amt) if not fp_closed else b""
    568     if amt is None:
    569         flush_decoder = True

File ~/code/climatematch/course-content/venv/lib/python3.8/site-packages/urllib3/response.py:533, in HTTPResponse._fp_read(self, amt)
    530     return buffer.getvalue()
    531 else:
    532     # StringIO doesn't like amt=None
--> 533     return self._fp.read(amt) if amt is not None else self._fp.read()

File /usr/lib/python3.8/http/client.py:459, in HTTPResponse.read(self, amt)
    456 if amt is not None:
    457     # Amount is given, implement using readinto
    458     b = bytearray(amt)
--> 459     n = self.readinto(b)
    460     return memoryview(b)[:n].tobytes()
    461 else:
    462     # Amount is not given (unbounded read) so we must check self.length
    463     # and self.chunked

File /usr/lib/python3.8/http/client.py:503, in HTTPResponse.readinto(self, b)
    498         b = memoryview(b)[0:self.length]
    500 # we do not use _safe_read() here because this may be a .will_close
    501 # connection, and the user is reading more bytes than will be provided
    502 # (for example, reading in 1k chunks)
--> 503 n = self.fp.readinto(b)
    504 if not n and b:
    505     # Ideally, we would raise IncompleteRead if the content-length
    506     # wasn't satisfied, but it might break compatibility.
    507     self._close_conn()

File /usr/lib/python3.8/socket.py:669, in SocketIO.readinto(self, b)
    667 while True:
    668     try:
--> 669         return self._sock.recv_into(b)
    670     except timeout:
    671         self._timeout_occurred = True

File /usr/lib/python3.8/ssl.py:1241, in SSLSocket.recv_into(self, buffer, nbytes, flags)
   1237     if flags != 0:
   1238         raise ValueError(
   1239           "non-zero flags not allowed in calls to recv_into() on %s" %
   1240           self.__class__)
-> 1241     return self.read(nbytes, buffer)
   1242 else:
   1243     return super().recv_into(buffer, nbytes, flags)

File /usr/lib/python3.8/ssl.py:1099, in SSLSocket.read(self, len, buffer)
   1097 try:
   1098     if buffer is not None:
-> 1099         return self._sslobj.read(len, buffer)
   1100     else:
   1101         return self._sslobj.read(len)

KeyboardInterrupt: 

ERA5-Land monthly averaged data from 1950 to present¶

ERA5-Land is a reanalysis dataset that offers an enhanced resolution compared to ERA5, providing a consistent view of land variables over several decades. It is created by replaying the land component of the ECMWF ERA5 climate reanalysis, which combines model data and global observations to generate a complete and reliable dataset using the laws of physics.

ERA5-Land focuses on the water and energy cycles at the surface level, offering a detailed record starting from 1950. The data used here is a post-processed subset of the complete ERA5-Land dataset. Monthly-mean averages have been pre-calculated to facilitate quick and convenient access to the data, particularly for applications that do not require sub-monthly fields. The native spatial resolution of the ERA5-Land reanalysis dataset is 9km on a reduced Gaussian grid (TCo1279). The data in the CDS has been regridded to a regular lat-lon grid of 0.1x0.1 degrees.

To calculate albedo using ERA5-Land¶

ERA5 parameter Forecast albedo provides is the measure of the reflectivity of the Earth’s surface. It is the fraction of solar (shortwave) radiation reflected by Earth’s surface, across the solar spectrum, for both direct and diffuse radiation. Values are between 0 and 1. Typically, snow and ice have high reflectivity with albedo values of 0.8 and above, land has intermediate values between about 0.1 and 0.4 and the ocean has low values of 0.1 or less. Radiation from the Sun (solar, or shortwave, radiation) is partly reflected back to space by clouds and particles in the atmosphere (aerosols) and some of it is absorbed. The rest is incident on the Earth’s surface, where some of it is reflected. The portion that is reflected by the Earth’s surface depends on the albedo. In the ECMWF Integrated Forecasting System (IFS), a climatological background albedo (observed values averaged over a period of several years) is used, modified by the model over water, ice and snow. Albedo is often shown as a percentage (%).

Climatematch Academy: Computational Tools for Climate Science

Change in Earth’s albedo and its dependence on land cover changes in the past 20 years

Contents