Mapping high-resolution Land Surface Temperature for UK cities¶
Description & purpose¶
This tutorial demonstrates how to access EOCIS Land Surface Temperature (LST) from Landsat 8 data from the Earth Observation Data Hub (EODH) using the pyeodh Python client. It covers catalogue searching, loading in cloud-optimised assets, and performing reprojection for mapping the data over national and city scales to produce interactive and static visualisations.
Setup 📦¶
Ensure the required Python libraries are available.
This includes the pyeodh client for interacting with the EODH data catalogue.
In [17]:
Copied!
# Run this cell if pyeodh is not installed, or needs updating
get_ipython().system('pip install --upgrade pyeodh')
get_ipython().system('pip install folium shapely xarray fsspec rioxarray cartopy pyproj')
import numpy as np
import pyeodh
import cartopy.crs as ccrs
import folium
import fsspec
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import rioxarray
import shapely as sh
import xarray as xr
from pyproj import Transformer
from shapely.geometry import Point
# Run this cell if pyeodh is not installed, or needs updating
get_ipython().system('pip install --upgrade pyeodh')
get_ipython().system('pip install folium shapely xarray fsspec rioxarray cartopy pyproj')
import numpy as np
import pyeodh
import cartopy.crs as ccrs
import folium
import fsspec
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import rioxarray
import shapely as sh
import xarray as xr
from pyproj import Transformer
from shapely.geometry import Point
Requirement already satisfied: pyeodh in /opt/conda/lib/python3.13/site-packages (0.1.12) Requirement already satisfied: lxml>=6.0.2 in /opt/conda/lib/python3.13/site-packages (from pyeodh) (6.0.2) Requirement already satisfied: owslib>=0.35.0 in /opt/conda/lib/python3.13/site-packages (from pyeodh) (0.35.0) Requirement already satisfied: pystac>=1.14.3 in /opt/conda/lib/python3.13/site-packages (from pyeodh) (1.14.3) Requirement already satisfied: requests-oauthlib>=2.0.0 in /opt/conda/lib/python3.13/site-packages (from pyeodh) (2.0.0) Requirement already satisfied: python-dateutil in /opt/conda/lib/python3.13/site-packages (from owslib>=0.35.0->pyeodh) (2.9.0.post0) Requirement already satisfied: pyyaml in /opt/conda/lib/python3.13/site-packages (from owslib>=0.35.0->pyeodh) (6.0.3) Requirement already satisfied: requests in /opt/conda/lib/python3.13/site-packages (from owslib>=0.35.0->pyeodh) (2.32.5) Requirement already satisfied: six>=1.5 in /opt/conda/lib/python3.13/site-packages (from python-dateutil->owslib>=0.35.0->pyeodh) (1.17.0) Requirement already satisfied: oauthlib>=3.0.0 in /opt/conda/lib/python3.13/site-packages (from requests-oauthlib>=2.0.0->pyeodh) (3.3.1) Requirement already satisfied: charset_normalizer<4,>=2 in /opt/conda/lib/python3.13/site-packages (from requests->owslib>=0.35.0->pyeodh) (3.4.4) Requirement already satisfied: idna<4,>=2.5 in /opt/conda/lib/python3.13/site-packages (from requests->owslib>=0.35.0->pyeodh) (3.11) Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/conda/lib/python3.13/site-packages (from requests->owslib>=0.35.0->pyeodh) (1.26.20) Requirement already satisfied: certifi>=2017.4.17 in /opt/conda/lib/python3.13/site-packages (from requests->owslib>=0.35.0->pyeodh) (2026.1.4) Requirement already satisfied: folium in /opt/conda/lib/python3.13/site-packages (0.20.0) Requirement already satisfied: shapely in /opt/conda/lib/python3.13/site-packages (2.1.0) Requirement already satisfied: xarray in /opt/conda/lib/python3.13/site-packages (2024.11.0) Requirement already satisfied: fsspec in /opt/conda/lib/python3.13/site-packages (2024.12.0) Requirement already satisfied: rioxarray in /opt/conda/lib/python3.13/site-packages (0.21.0) Requirement already satisfied: cartopy in /opt/conda/lib/python3.13/site-packages (0.25.0) Requirement already satisfied: pyproj in /opt/conda/lib/python3.13/site-packages (3.7.2) Requirement already satisfied: branca>=0.6.0 in /opt/conda/lib/python3.13/site-packages (from folium) (0.8.2) Requirement already satisfied: jinja2>=2.9 in /opt/conda/lib/python3.13/site-packages (from folium) (3.1.6) Requirement already satisfied: numpy in /opt/conda/lib/python3.13/site-packages (from folium) (2.3.5) Requirement already satisfied: requests in /opt/conda/lib/python3.13/site-packages (from folium) (2.32.5) Requirement already satisfied: xyzservices in /opt/conda/lib/python3.13/site-packages (from folium) (2025.11.0) Requirement already satisfied: packaging>=23.2 in /opt/conda/lib/python3.13/site-packages (from xarray) (26.0) Requirement already satisfied: pandas>=2.1 in /opt/conda/lib/python3.13/site-packages (from xarray) (3.0.1) Requirement already satisfied: rasterio>=1.4.3 in /opt/conda/lib/python3.13/site-packages (from rioxarray) (1.5.0) Requirement already satisfied: matplotlib>=3.6 in /opt/conda/lib/python3.13/site-packages (from cartopy) (3.10.8) Requirement already satisfied: pyshp>=2.3 in /opt/conda/lib/python3.13/site-packages (from cartopy) (3.1.6) Requirement already satisfied: certifi in /opt/conda/lib/python3.13/site-packages (from pyproj) (2026.1.4) Requirement already satisfied: MarkupSafe>=2.0 in /opt/conda/lib/python3.13/site-packages (from jinja2>=2.9->folium) (3.0.3) Requirement already satisfied: contourpy>=1.0.1 in /opt/conda/lib/python3.13/site-packages (from matplotlib>=3.6->cartopy) (1.3.3) Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.13/site-packages (from matplotlib>=3.6->cartopy) (0.12.1) Requirement already satisfied: fonttools>=4.22.0 in /opt/conda/lib/python3.13/site-packages (from matplotlib>=3.6->cartopy) (4.61.1) Requirement already satisfied: kiwisolver>=1.3.1 in /opt/conda/lib/python3.13/site-packages (from matplotlib>=3.6->cartopy) (1.4.9) Requirement already satisfied: pillow>=8 in /opt/conda/lib/python3.13/site-packages (from matplotlib>=3.6->cartopy) (12.1.1) Requirement already satisfied: pyparsing>=3 in /opt/conda/lib/python3.13/site-packages (from matplotlib>=3.6->cartopy) (3.3.2) Requirement already satisfied: python-dateutil>=2.7 in /opt/conda/lib/python3.13/site-packages (from matplotlib>=3.6->cartopy) (2.9.0.post0) Requirement already satisfied: six>=1.5 in /opt/conda/lib/python3.13/site-packages (from python-dateutil>=2.7->matplotlib>=3.6->cartopy) (1.17.0) Requirement already satisfied: affine in /opt/conda/lib/python3.13/site-packages (from rasterio>=1.4.3->rioxarray) (2.4.0) Requirement already satisfied: attrs in /opt/conda/lib/python3.13/site-packages (from rasterio>=1.4.3->rioxarray) (25.4.0) Requirement already satisfied: click!=8.2.*,>=4.0 in /opt/conda/lib/python3.13/site-packages (from rasterio>=1.4.3->rioxarray) (8.3.1) Requirement already satisfied: cligj>=0.5 in /opt/conda/lib/python3.13/site-packages (from rasterio>=1.4.3->rioxarray) (0.7.2) Requirement already satisfied: charset_normalizer<4,>=2 in /opt/conda/lib/python3.13/site-packages (from requests->folium) (3.4.4) Requirement already satisfied: idna<4,>=2.5 in /opt/conda/lib/python3.13/site-packages (from requests->folium) (3.11) Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/conda/lib/python3.13/site-packages (from requests->folium) (1.26.20)
Connect to the EODH 🔗¶
We first create an instance of the Python client and connect to the catalogue service. A time range and spatial extent are also defined for use throughout the analysis.
In [6]:
Copied!
# Connect to the Hub
client = pyeodh.Client(base_url="https://eodatahub.org.uk").get_catalog_service()
# Define temporal and spatial range for the collection search.
# Selected the 26th of August 2019 as it was a clear day according to historical MET data.
time_start = "2019-08-26"
time_end = "2019-08-26"
#london_slice
lat_range = [51.3100, 51.6700]
lon_range = [-0.5200, 0.2100]
# Connect to the Hub
client = pyeodh.Client(base_url="https://eodatahub.org.uk").get_catalog_service()
# Define temporal and spatial range for the collection search.
# Selected the 26th of August 2019 as it was a clear day according to historical MET data.
time_start = "2019-08-26"
time_end = "2019-08-26"
#london_slice
lat_range = [51.3100, 51.6700]
lon_range = [-0.5200, 0.2100]
Search the catalogue 🔍¶
The EODH catalogue can be searched programmatically using dataset identifiers and metadata filters.
Here we search for Daily Sentinel‑3B Land Surface Temperature products within the selected time window, using the eocis-chuk-lst identifier.
In [7]:
Copied!
datasearch = client.search(
collections=["eocis-chuk-lst"],
catalog_paths=["public/catalogs/ceda-stac-catalogue"],
start_datetime=time_start,
end_datetime=time_end,
)
datasearch = client.search(
collections=["eocis-chuk-lst"],
catalog_paths=["public/catalogs/ceda-stac-catalogue"],
start_datetime=time_start,
end_datetime=time_end,
)
Each search result may contain multiple assets. Where available, Kerchunk reference files are preferred, as they allow efficient remote access to NetCDF datasets without downloading the full data volume.
In [10]:
Copied!
list_of_data = []
for item in datasearch:
asset = None
if "reference_file" in item.assets:
asset = item.assets["reference_file"]
href = asset.href
is_kerchunk = True
else:
for k, a in item.assets.items():
if "data" in a.roles or a.type == "application/netcdf":
asset = a
href = asset.href
is_kerchunk = False
break
if is_kerchunk:
list_of_data.append(href)
else:
print("nc")
list_of_data = []
for item in datasearch:
asset = None
if "reference_file" in item.assets:
asset = item.assets["reference_file"]
href = asset.href
is_kerchunk = True
else:
for k, a in item.assets.items():
if "data" in a.roles or a.type == "application/netcdf":
asset = a
href = asset.href
is_kerchunk = False
break
if is_kerchunk:
list_of_data.append(href)
else:
print("nc")
Load in the data 🌐¶
The collected Kerchunk reference is opened using xarray enabling time‑aware and spatial operations.
In [12]:
Copied!
ds = xr.open_dataset(list_of_data[0], engine="kerchunk")
ds
ds = xr.open_dataset(list_of_data[0], engine="kerchunk")
ds
Out[12]:
<xarray.Dataset> Size: 666MB
Dimensions: (y: 15170, x: 10970, bnds: 2)
Coordinates:
* x (x) int32 44kB -331950 -331850 -331750 ... 764750 764850 764950
* y (y) int32 61kB 1249950 1249850 1249750 ... -266750 -266850 -266950
Dimensions without coordinates: bnds
Data variables:
crsOSGB float64 8B ...
lst (y, x) float32 666MB ...
x_bnds (x, bnds) int32 88kB ...
y_bnds (y, bnds) int32 121kB ...
Attributes: (12/41)
Convention: CF-1.8
Conventions: CF-1.10
acknowledgement: Funded by the Natural Environment Research Co...
cdm_data_type: grid
comment: This dataset includes regridded data from Lan...
creator_email:
... ...
time_coverage_duration: P1D
time_coverage_end: 20190826T235959Z
time_coverage_resolution: P1D
time_coverage_start: 20190826T000000Z
title: EOCIS CHUK Daily Daytime Land Surface Tempera...
version: 1.0Here we slice the dataset by subsetting the array to a UK bounding box, using latitude and longitude coordinates we defined earlier. This reduces the amount of data we are dealing with, storing only what we need for the visualisation.
Since this dataset is stored in BNG x and y points, we run a quick loop-up to match the specified lat/lons to the indexes in the file.
In [13]:
Copied!
transformer = Transformer.from_crs(
"EPSG:4326", # lon/lat
"EPSG:27700", # British National Grid
always_xy=True
)
# Convert all four corners
corners = [
(lon_range[0], lat_range[0]), # SW
(lon_range[1], lat_range[0]), # SE
(lon_range[0], lat_range[1]), # NW
(lon_range[1], lat_range[1]), # NE
]
xy = [transformer.transform(lon, lat) for lon, lat in corners]
eastings = [x for x, y in xy]
northings = [y for x, y in xy]
x_min, x_max = min(eastings), max(eastings)
y_min, y_max = min(northings), max(northings)
cutout = ds.sel(
x=slice(x_min, x_max),
y=slice(y_max, y_min)
)
transformer = Transformer.from_crs(
"EPSG:4326", # lon/lat
"EPSG:27700", # British National Grid
always_xy=True
)
# Convert all four corners
corners = [
(lon_range[0], lat_range[0]), # SW
(lon_range[1], lat_range[0]), # SE
(lon_range[0], lat_range[1]), # NW
(lon_range[1], lat_range[1]), # NE
]
xy = [transformer.transform(lon, lat) for lon, lat in corners]
eastings = [x for x, y in xy]
northings = [y for x, y in xy]
x_min, x_max = min(eastings), max(eastings)
y_min, y_max = min(northings), max(northings)
cutout = ds.sel(
x=slice(x_min, x_max),
y=slice(y_max, y_min)
)
In [11]:
Copied!
cutout.lst.plot()
cutout.lst.plot()
Out[11]:
<matplotlib.collections.QuadMesh at 0x7f691814c2f0>
In [18]:
Copied!
X, Y = np.meshgrid(cutout.x.values,cutout.y.values)
transformer = Transformer.from_crs(
"EPSG:27700",
"EPSG:4326",
always_xy=True
)
lon, lat = transformer.transform(X, Y)
X, Y = np.meshgrid(cutout.x.values,cutout.y.values)
transformer = Transformer.from_crs(
"EPSG:27700",
"EPSG:4326",
always_xy=True
)
lon, lat = transformer.transform(X, Y)
In [45]:
Copied!
proj = ccrs.OSGB()
plt.figure(figsize=(10,10))
ax = plt.axes(projection=proj)
ax.gridlines(draw_labels=True,zorder=0)
mesh = cutout.lst.plot(
ax=ax,
transform=proj, # IMPORTANT: matches data CRS
cmap="inferno",
add_colorbar=False,
#cbar_kwargs={"label": "Land Surface Temperature (°C)"}
zorder=3
)
ax_ins = ax.inset_axes([0.8, 0.15, 0.95, 0.2], projection=proj)
ax_ins.axis('off')
cb = plt.colorbar(mesh, ax=ax_ins, location="left", pad=0, aspect=10)
cb.ax.set_title("LST [K]")
# Set the background color and alpha
cb.ax.set_facecolor('lightgrey')
cb.ax.patch.set_alpha(0.5) # Make it semi-transparent so colormap shows through
# Move the patch forward so it sits on top of the colormap mesh
cb.ax.patch.set_zorder(-1)
cb.ax.collections[0].set_zorder(0)
#ax_ins.colorbar(mesh, orientation='vertical')
#ax.coastlines()
ax.set_title("Land Surface Temperature 2019-08-26")
plt.show()
proj = ccrs.OSGB()
plt.figure(figsize=(10,10))
ax = plt.axes(projection=proj)
ax.gridlines(draw_labels=True,zorder=0)
mesh = cutout.lst.plot(
ax=ax,
transform=proj, # IMPORTANT: matches data CRS
cmap="inferno",
add_colorbar=False,
#cbar_kwargs={"label": "Land Surface Temperature (°C)"}
zorder=3
)
ax_ins = ax.inset_axes([0.8, 0.15, 0.95, 0.2], projection=proj)
ax_ins.axis('off')
cb = plt.colorbar(mesh, ax=ax_ins, location="left", pad=0, aspect=10)
cb.ax.set_title("LST [K]")
# Set the background color and alpha
cb.ax.set_facecolor('lightgrey')
cb.ax.patch.set_alpha(0.5) # Make it semi-transparent so colormap shows through
# Move the patch forward so it sits on top of the colormap mesh
cb.ax.patch.set_zorder(-1)
cb.ax.collections[0].set_zorder(0)
#ax_ins.colorbar(mesh, orientation='vertical')
#ax.coastlines()
ax.set_title("Land Surface Temperature 2019-08-26")
plt.show()
In [47]:
Copied!
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import cartopy.crs as ccrs
proj = ccrs.OSGB()
fig = plt.figure(figsize=(10,10))
ax = plt.axes(projection=proj)
ax.gridlines(draw_labels=True, zorder=0)
# Plot data
mesh = cutout.lst.plot(
ax=ax,
transform=proj,
cmap="inferno",
add_colorbar=False,
zorder=3
)
# 1. Create a standard 2D inset axis to serve as our solid background block.
# [left, bottom, width, height] relative to main plot.
cax_bg = ax.inset_axes([0.75, 0.15, 0.20, 0.25])
cax_bg.set_facecolor('lightgrey') # This will be your solid background fill
cax_bg.set_xticks([]) # Hide background axes ticks
cax_bg.set_yticks([])
# 2. Nest a smaller, tight colorbar axis INSIDE that background block.
# [left, bottom, width, height] relative to the cax_bg frame itself.
# We leave room on the left (0.5) for the labels/ticks, and room at the top (0.85) for the title.
cax = cax_bg.inset_axes([0.6, 0.08, 0.15, 0.75])
# 3. Direct the colorbar to use our explicit custom axis via 'cax'
cb = plt.colorbar(mesh, cax=cax, orientation='vertical')
# 4. Move ticks to the left side and format text
cax.yaxis.set_ticks_position('left')
cb.ax.set_title("LST [K]", pad=10, fontsize=10)
ax.set_title("Land Surface Temperature 2019-08-26")
plt.show()
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import cartopy.crs as ccrs
proj = ccrs.OSGB()
fig = plt.figure(figsize=(10,10))
ax = plt.axes(projection=proj)
ax.gridlines(draw_labels=True, zorder=0)
# Plot data
mesh = cutout.lst.plot(
ax=ax,
transform=proj,
cmap="inferno",
add_colorbar=False,
zorder=3
)
# 1. Create a standard 2D inset axis to serve as our solid background block.
# [left, bottom, width, height] relative to main plot.
cax_bg = ax.inset_axes([0.75, 0.15, 0.20, 0.25])
cax_bg.set_facecolor('lightgrey') # This will be your solid background fill
cax_bg.set_xticks([]) # Hide background axes ticks
cax_bg.set_yticks([])
# 2. Nest a smaller, tight colorbar axis INSIDE that background block.
# [left, bottom, width, height] relative to the cax_bg frame itself.
# We leave room on the left (0.5) for the labels/ticks, and room at the top (0.85) for the title.
cax = cax_bg.inset_axes([0.6, 0.08, 0.15, 0.75])
# 3. Direct the colorbar to use our explicit custom axis via 'cax'
cb = plt.colorbar(mesh, cax=cax, orientation='vertical')
# 4. Move ticks to the left side and format text
cax.yaxis.set_ticks_position('left')
cb.ax.set_title("LST [K]", pad=10, fontsize=10)
ax.set_title("Land Surface Temperature 2019-08-26")
plt.show()
In [ ]:
Copied!