Reading GHCN observations

This notebook introduces how to use the hydropandas package to read, process and visualise data from the Global Historical Climatology Network (GHCN). GHCN-Daily is a database of daily climate summaries from land surface stations around the world, maintained by NOAA.

[1]:

import contextily as ctx
import matplotlib.pyplot as plt

import hydropandas as hpd

# enabling logging so we can see what happens in the background
hpd.util.get_color_logger("INFO");

Read GHCN observations within an extent

Use hpd.read_ghcn to download all GHCN stations within a bounding box. The extent parameter is [xmin, xmax, ymin, ymax] in the coordinate system given by crs. Setting elements='PRCP' downloads only precipitation data; leave it as None to retrieve all available elements.

[2]:

# read GHCN precipitation observations for an area in the Netherlands
# extent: [xmin, xmax, ymin, ymax] in RD New (EPSG:28992)
extent = [130_000, 150_000, 450_000, 470_000]

oc = hpd.read_ghcn(
    extent=extent,
    crs=28992,
    elements="PRCP",
    tmin="2020-01-01",
    tmax="2021-12-31",
)
oc

INFO:hydropandas.io.ghcn.get_obs_list_from_extent:downloading GHCN data from 3 stations

station: 100%|██████████| 3/3 [00:02<00:00,  1.47it/s]

[2]:

	meteo_var	x	source	unit	y	station	obs
name
NLE00101997_PRCP	PRCP	130800.683837	GHCN	m	469904.652522	NLE00101997	MeteoObs NLE00101997_PRCP -----metadata------ ...
NLE00152487		147199.951428	GHCN		459715.887745	NLE00152487	MeteoObs NLE00152487 -----metadata------ meteo...
NLM00006260_PRCP	PRCP	140761.236278	GHCN	m	456759.504225	NLM00006260	MeteoObs NLM00006260_PRCP -----metadata------ ...

[3]:

# plot station locations on a map
gdf = oc.to_gdf().to_crs(epsg=4326)
ax = gdf.plot(figsize=(8, 8), color="steelblue", markersize=60)
ctx.add_basemap(ax=ax, crs=4326, attribution=False)

for idx, row in gdf.iterrows():
    ax.annotate(
        text=idx,
        xy=(row.geometry.x, row.geometry.y),
        fontsize=7,
        ha="center",
        va="bottom",
    )
ax.set_title("GHCN stations")
plt.tight_layout()

Plot observations from a single station

[4]:

# select the first observation and plot the precipitation time series
o = oc.iloc[0].obs
print(f"Station: {o.name}  |  element: PRCP  |  unit: {o.unit}")
o["PRCP"].plot(
    figsize=(12, 4),
    drawstyle="steps",
    ylabel="Precipitation (m)",
    title=f"Daily precipitation – {o.name}",
)
plt.tight_layout()

Station: NLE00101997_PRCP  |  element: PRCP  |  unit: m

Download multiple elements (temperature)

Omit the elements argument (or pass a list) to download all available elements for the stations in the extent. Here we request maximum and minimum daily temperature (TMAX and TMIN) for a small extent around De Bilt.

[5]:

# read TMAX and TMIN for De Bilt area
extent_debilt = [4.9, 5.2, 51.9, 52.1]  # [xmin, xmax, ymin, ymax] WGS84

oc_temp = hpd.read_ghcn(
    extent=extent_debilt,
    crs=4326,
    elements=["TMAX", "TMIN"],
    tmin="2020-01-01",
    tmax="2020-12-31",
)
oc_temp

INFO:hydropandas.io.ghcn.get_obs_list_from_extent:downloading GHCN data from 4 stations

station: 100%|██████████| 4/4 [00:02<00:00,  1.60it/s]

[5]:

	meteo_var	x	source	unit	y	station	obs
name
NLE00108992		4.9331	GHCN		52.0000	NLE00108992	MeteoObs NLE00108992 -----metadata------ meteo...
NLE00109338		5.0000	GHCN		52.0831	NLE00109338	MeteoObs NLE00109338 -----metadata------ meteo...
NLE00152476_TMAX	TMAX	4.9256	GHCN	unknown	51.9692	NLE00152476	MeteoObs NLE00152476_TMAX -----metadata------ ...
NLE00152476_TMIN	TMIN	4.9256	GHCN	unknown	51.9692	NLE00152476	MeteoObs NLE00152476_TMIN -----metadata------ ...
NLM00006260_TMAX	TMAX	5.1794	GHCN	unknown	52.0989	NLM00006260	MeteoObs NLM00006260_TMAX -----metadata------ ...
NLM00006260_TMIN	TMIN	5.1794	GHCN	unknown	52.0989	NLM00006260	MeteoObs NLM00006260_TMIN -----metadata------ ...

[6]:

# plot temperature observations for the first station
temp_min = oc_temp.get_obs(meteo_var="TMIN", station="NLE00152476")
temp_max = oc_temp.get_obs(meteo_var="TMAX", station="NLE00152476")
f, ax = plt.subplots(figsize=(12, 4))
temp_max["TMAX"].plot(ax=ax, label="TMAX", color="tomato")
temp_min["TMIN"].plot(ax=ax, label="TMIN", color="steelblue")
ax.set_ylabel("Temperature (°C / 0.1 °C)")
ax.set_title(f"Daily temperature – {temp_min.station}")
ax.legend()
f.tight_layout()

Explore available stations

Use ghcn.get_stations to retrieve a GeoDataFrame with all GHCN stations worldwide or filtered to an extent, without downloading any measurements.

[7]:

# get station metadata only (fast – no measurements downloaded)
oc_meta = hpd.read_ghcn(
    extent=extent_debilt,
    crs=4326,
    only_metadata=True,
)
print(f"Found {len(oc_meta)} stations in the extent")
oc_meta[["x", "y", "source", "station"]]

INFO:hydropandas.io.ghcn.get_obs_list_from_extent:downloading GHCN data from 4 stations

station: 100%|██████████| 4/4 [00:00<00:00, 1982.42it/s]

Found 4 stations in the extent

[7]:

	x	y	source	station
name
NLE00108992	4.9331	52.0000	GHCN	NLE00108992
NLE00109338	5.0000	52.0831	GHCN	NLE00109338
NLE00152476	4.9256	51.9692	GHCN	NLE00152476
NLM00006260	5.1794	52.0989	GHCN	NLM00006260

[8]:

# interactive map of all downloaded observations
oc.plots.interactive_map(plot_dir="figure", per_location=False, popup_width=300)

INFO:hydropandas.extensions.plots.interactive_map:no iplot available for NLE00152487

[8]:

Make this Notebook Trusted to load map: File -> Trust Notebook