calculate_heat_hazard_indicators.py

"""
Generate indicator variables for South Gloucestershire heat hazard analysis

James Thomas, Eunice Lo
March 2023

Using dataset: https://catalogue.ceda.ac.uk/uuid/bbca3267dc7d4219af484976734c9527
Citable as: Met Office; Hollis, D.; McCarthy, M.; Kendon, M.; Legg, T. (2022): HadUK-Grid Gridded Climate Observations on a 1km grid over the UK, v1.1.0.0 (1836-2021). NERC EDS Centre for Environmental Data Analysis, 26 May 2022. doi:10.5285/bbca3267dc7d4219af484976734c9527. https://dx.doi.org/10.5285/bbca3267dc7d4219af484976734c9527

Processed on JASMIN using data from the CEDA Archive (in approx 1 minute)
Acknowledge as: This work was carried out using the computational facilities of JASMIN – https://www.jasmin.ac.uk/, using data hosted by the CEDA Archive - https://archive.ceda.ac.uk/.

Also processed on UoBristol HPC
Acknowledge as: This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol - http://www.bristol.ac.uk/acrc/.
"""

import glob

import numpy as np
from pyproj import Transformer
import xarray as xr


# Locate the datafiles of interest:
# Summer months (JJA) for 2015-2019

datafile_template = '/badc/ukmo-hadobs/data/insitu/MOHC/HadOBS/HadUK-Grid/v1.1.0.0/1km/{variable}/day/v20220310/{variable}_hadukgrid_uk_1km_day_{year}{month:02d}01-*.nc'

datafile_paths = []
for variable in ['tasmax', 'tasmin']:
    for year in range(2015, 2020):
        for month in range(6, 9):
            datafile_pattern = datafile_template.format(variable=variable, year=year, month=month)
            datafile_paths.extend(glob.glob(datafile_pattern))


# Open the datafiles and select on the area of interest

dataset = xr.open_mfdataset(datafile_paths).rio.set_crs(27700)

wgs84_to_osgb36 = Transformer.from_crs(4326, 27700, always_xy=True)          # Care with argument order
left, bottom, right, top = wgs84_to_osgb36.transform_bounds(-3, 51, -2, 52)  # (using always_xy=True)
study_area = dataset.sel(projection_x_coordinate=slice(left, right), projection_y_coordinate=slice(bottom, top))


# Calculate the statistics of interest

# Average daily maximum temperature
average_daily_maximum = (
    study_area
    .tasmax
    .mean(dim='time')
    .rename('average_daily_maximum_temperature')
)

# Proportion of days where "average" temperature above a threshold
# (where average is estimated from the mean of min and max daily temperature)
tas = (study_area.tasmax + study_area.tasmin) / 2
threshold = 15
days_over_threshold = (
    (tas > 15)
    .mean(dim='time')
    .where(~average_daily_maximum.isnull())  # Mask areas with no data
    .rename('days_average_temperature_over_threshold')
)


# Output data in netCDF and GeoTIFF formats

# Re-set the CRS because it gets forgotten
average_daily_maximum = average_daily_maximum.rio.set_crs(27700)
days_over_threshold = days_over_threshold.rio.set_crs(27700)

average_daily_maximum.to_netcdf('average_daily_maximum.nc')
average_daily_maximum.rio.to_raster('average_daily_maximum.tif')
days_over_threshold.to_netcdf('days_over_threshold.nc')
days_over_threshold.rio.to_raster('days_over_threshold.tif')