Two numbers that people often want for drought, heat and crop work are how warm the ground is and how wet the soil is. They are not the same as air temperature and rainfall. Land surface temperature, LST, is how hot the ground surface is, as a satellite sees it, and on a summer afternoon it can be 15 degrees hotter than the air above it. Soil moisture is the share of the soil volume that is water, in cubic metres of water per cubic metre of soil, so a value of 0.20 means that a fifth of the volume is water.
The two come from different places. LST is measured by satellites, and this guide reads the MODIS product, at 1 km. Soil moisture on a long record, for any place, comes from a weather model that mixes a land model with observations, and this guide reads it from the Open-Meteo archive, which uses ERA5-Land. Satellite soil moisture products such as NASA SMAP need a free NASA Earthdata account, and they are not used here.
This guide has one route, in Python, because both sources are web services with no account. The study place is Ludhiana, Punjab, and the code uses ludhiana_district.geojson from the shared data of this section.
Soil moisture for a place
You need Python with requests, pandas, numpy, rasterio and geopandas.
import pandas as pd
import requests
# Step 1: ask Open-Meteo for the hourly soil moisture of the top 7 cm and the next layer, and the rain, for 2023. Ludhiana city.
params = {
"latitude": 30.90, "longitude": 75.85,
"start_date": "2023-01-01", "end_date": "2023-12-31",
"hourly": "soil_moisture_0_to_7cm,soil_moisture_7_to_28cm,precipitation",
"timezone": "Asia/Kolkata",
}
reply = requests.get("https://archive-api.open-meteo.com/v1/archive", params=params, headers={"User-Agent": "timeandspatial-guide/1.0 (write your email address here)"}, timeout=120)
reply.raise_for_status()
reply.encoding = "utf-8"
data = reply.json()
print("grid cell used:", round(data["latitude"], 2), round(data["longitude"], 2), "| units:", data["hourly_units"]["soil_moisture_0_to_7cm"].replace("³", "3"))
print("hours returned:", len(data["hourly"]["time"]))
grid cell used: 30.9 75.83 | units: m3/m3
hours returned: 8760
# Step 1: put the hourly numbers in a table with a time index
hourly = pd.DataFrame(data["hourly"])
hourly["time"] = pd.to_datetime(hourly["time"])
hourly = hourly.set_index("time")
# Step 2: monthly means for the soil moisture and monthly totals for the rain
monthly = pd.DataFrame({
"soil 0-7 cm": hourly["soil_moisture_0_to_7cm"].resample("MS").mean(),
"soil 7-28 cm": hourly["soil_moisture_7_to_28cm"].resample("MS").mean(),
"rain mm": hourly["precipitation"].resample("MS").sum(),
})
monthly.index = monthly.index.strftime("%b")
print(monthly.round(3).T.to_string())
print("rain for the year:", round(float(monthly["rain mm"].sum())), "mm")
time Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
soil 0-7 cm 0.292 0.220 0.268 0.221 0.217 0.240 0.370 0.325 0.245 0.206 0.177 0.183
soil 7-28 cm 0.290 0.234 0.250 0.245 0.208 0.243 0.357 0.316 0.241 0.227 0.201 0.193
rain mm 49.000 4.100 81.000 25.800 113.100 51.600 211.900 92.800 80.400 18.500 11.700 0.000
rain for the year: 740 mm
The hourly values are in cubic metres per cubic metre, written m3/m3 above. In the table, the top layer is wettest in July, at 0.37, the month with the most rain, 212 mm, and driest in November, at 0.18. The top layer follows the rain and drops quickly when the rain stops, because the top 7 cm dries in days. The layer below, 7 to 28 cm, moves more slowly, and in the dry months of October, November and December it is wetter than the top layer, by 0.01 to 0.02. The rain in this table comes from the same model, ERA5-Land, so it is not the same as the CHIRPS rain of the rainfall guide.
The grid cell that answers is about 9 km wide, as Open-Meteo says in its documentation. The number describes that cell, and not a field.
Land surface temperature from MODIS
MODIS flies on two satellites. Terra passes over Punjab at about 10:30 in the morning and Aqua at about 13:30 in the afternoon, and the daily products are called MOD11A1 for Terra and MYD11A1 for Aqua. A pixel is 1 km. The product stores the temperature in kelvin multiplied by 50, so you multiply by 0.02 and subtract 273.15 to get degrees Celsius. A value of 0 means no data, usually because of cloud.
The Planetary Computer keeps the MODIS files as cloud optimised GeoTIFFs and lets you read them with a short-lived token that it gives away without an account.
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.mask import mask
from rasterio.warp import transform_geom
# Step 1: the district, and the catalogue address
district = gpd.read_file("ludhiana_workflow_data/ludhiana_district.geojson")
STAC = "https://planetarycomputer.microsoft.com/api/stac/v1"
# Step 2: search the daily Aqua product (afternoon pass) for May and June 2023 over the district
search = {
"collections": ["modis-11A1-061"],
"bbox": [float(v) for v in district.total_bounds],
"datetime": "2023-05-01T00:00:00Z/2023-06-30T23:59:59Z",
"limit": 200,
}
items = requests.post(STAC + "/search", json=search, timeout=120).json()["features"]
aqua = sorted((i for i in items if i["id"].startswith("MYD11A1") and ".h24v05." in i["id"]), key=lambda i: i["properties"]["start_datetime"])
print(len(items), "items found,", len(aqua), "from Aqua in tile h24v05")
print("first:", aqua[0]["id"])
# Step 3: a token to read the files
token = requests.get("https://planetarycomputer.microsoft.com/api/sas/v1/token/modiseuwest/modis-061-cogs", timeout=60).json()["token"]
127 items found, 66 from Aqua in tile h24v05
first: MYD11A1.A2023121.h24v05.061.2023123210321
def district_lst(item, band="LST_Day_1km"):
# Step 1: open the band, cut it to the district, and convert to degrees Celsius
with rasterio.open(item["assets"][band]["href"] + "?" + token) as src:
shape = transform_geom("EPSG:4326", src.crs, district.geometry.iloc[0].__geo_interface__)
arr, _ = mask(src, [shape], crop=True, nodata=0)
arr = arr[0].astype("float64")
arr[arr == 0] = np.nan # 0 means no data
return arr * 0.02 - 273.15
# Step 2: the mean of the district on each Aqua day that has enough valid cells
rows = []
for item in aqua:
lst = district_lst(item)
valid = np.isfinite(lst)
if valid.sum() < 2000: # skip days when about half of the district or more is cloud
continue
rows.append({"date": item["properties"]["start_datetime"][:10], "valid cells": int(valid.sum()), "mean C": round(float(np.nanmean(lst)), 1), "max C": round(float(np.nanmax(lst)), 1)})
table = pd.DataFrame(rows).set_index("date")
print(len(table), "usable days of", len(aqua))
print(table.head(6).to_string())
# Step 3: the monthly mean of the usable days, next to the soil moisture and the rain of the same month
table.index = pd.to_datetime(table.index)
may_june = table["mean C"].resample("MS").mean().round(1)
may_june.index = may_june.index.strftime("%b")
side = pd.DataFrame({"LST mean C": may_june, "usable days": table["mean C"].resample("MS").count().to_numpy()})
side["soil 0-7 cm"] = monthly.loc[side.index, "soil 0-7 cm"].round(3)
side["rain mm"] = monthly.loc[side.index, "rain mm"].round(0)
print(side.to_string())
41 usable days of 66
valid cells mean C max C
date
2023-05-04 4293 41.6 47.9
2023-05-05 4258 37.1 40.4
2023-05-06 4277 42.7 48.1
2023-05-07 4316 37.9 43.4
2023-05-08 4274 36.1 40.2
2023-05-09 4317 45.7 50.7
LST mean C usable days soil 0-7 cm rain mm
date
May 42.3 27 0.217 113.0
Jun 42.0 14 0.240 52.0
MODIS gives the surface temperature of the ground, trees and roofs, not of the air. In the table, the district mean on a May afternoon runs from 36 to 46 degrees, and the hottest pixels reach 50.7 degrees. The last table puts the pieces side by side. Of the 66 Aqua days in May and June, 41 were usable, 27 in May and 14 in June. The mean afternoon surface temperature is 42.3 degrees in May and 42.0 degrees in June. The soil is a little wetter in June, 0.24 against 0.22, and the rain of the model is 113 mm in May and 52 mm in June. Two months are too few to say how the temperature follows the soil moisture, and for that you would run the same code for a whole year. A day with few valid cells is mostly cloud, and its mean is not a fair picture of the district, which is why the code drops days with fewer than 2000 valid cells, about half of the district.
What these numbers cannot tell you
- LST is the temperature of the surface, and it depends on the cover. Bare dry soil and a roof are hot, and a wet field or a canopy is cooler.
- A MODIS pixel is 1 km. A park or a canal inside a city is lost in it.
- MODIS only sees through clear sky. In the monsoon, many days have no usable value.
- The soil moisture of Open-Meteo is a model result for a 9 km cell. It is not a measurement of one field, and it does not know about irrigation, which in Punjab can make a field much wetter than the model says.
- Soil moisture of the top 7 cm changes within a day, so one hour tells you little.
Common mistakes
- Reading the MODIS number as degrees. It is kelvin times 50, and 0 means no data.
- Taking a mean that includes the 0 cells.
- Comparing a Terra day and an Aqua day. The overpass times differ by three hours, and so do the temperatures.
- Reading soil moisture as a percentage of rain. It is a share of the soil volume.
- Quoting the Open-Meteo rain as measured rain. It is a model value.