GB potential evapotranspiration dataset

A few weeks ago I wrote about a method to derive incoming solar radiation for different latitudes. This was to help a GIS Masters student who’s looking to derive an equation for propwet that can be used to show the effects of catchment management changes. As an intermediate step the student also needed to calculate potential evapotranspiration. He did this on a catchment basis for Scotland. PE data is commercially available in the UK from the Met Office as part of their MORECS data.

Feeling inspired from my radiation success I’ve extended my previous efforts to calculate PE on a 5 km grid for Great Britain. I followed the methodology outlined by Kay and Davies, which only required temperature data (section 2.1.2) – extract below. To make life easier I followed the suggestion that sensible results can be obtained on a monthly aggregation. For temperature I downloaded the 1969-1999 long term average data available from the Met Office as part of their UKCIP09 assessment. This comes in a raster format, one grid for each month.

Extract from Kay, A. L.; Davies, H. N.. 2008 Calculating potential evaporation from climate model data: a source of uncertainty for hydrological climate change impacts. Journal of Hydrology, 358 (3-4). 221-239.

Extract from Kay, A. L.; Davies, H. N.. 2008 Calculating potential evaporation from climate model data: a source of uncertainty for hydrological climate change impacts. Journal of Hydrology, 358 (3-4). 221-239.

The early grid and centroid steps taken are as described on my external radiation post, beyond that I carried out the following operations in QGIS to prepare the data:

  • Buffered an outline of the GB coast (vector/geoprocessing/buffer)
  • Selected centroids only within the buffer (vector/data management tools/join by location) to ease the processing time
  • Used the point sampling tool plugin to sample the 12 temperature raster layers and named the columns jan, feb, mar, etc.
  • Converted the British National Grid projection to a WGS84 (right click layer, save as)
  • Added lat/long data to the attribute table (vector/geometry tools/add geometry columns)
  • Tidied up the table to have a column order of x,y,jan,feb,mar… etc. (table manager plugin)
  • exported result as a csv file (right click layer, save as).

From this point I worked in R. I find the RStudio package very useful and do my code writing in RStudio or geany.

# import latitude
lat_raw = as.matrix(read.csv("~/filepath/MO_grid_pts.csv",sep=",",header=1))
# strip null temp data
lat_raw_na=na.omit(lat_raw)

# prepare latitude
lat_len = nrow(lat_raw_na)
lat_rad = lat_raw_na[ ,2]/180*pi
lat = matrix(lat_rad, nrow=lat_len, ncol=365)

# matrix of days, row number from lat
j = matrix(c(1:365), nrow=lat_len, ncol=365, byrow=1)

# calculate extraterrestial radiation in Mj/m2/day
sigma = .409*sin(.0172*j-1.39)
ws = acos(-tan(lat)*tan(sigma))
dr = 1+.033*cos(.0172*j)
Ra = 37.6*dr*(ws*sin(lat)*sin(sigma)+cos(lat)*cos(sigma)*sin(ws))

# mean Re by months
jan = rowMeans(Ra[,1:31])
feb = rowMeans(Ra[,32:59])
mar = rowMeans(Ra[,60:90])
apr = rowMeans(Ra[,91:120])
may = rowMeans(Ra[,121:151])
jun = rowMeans(Ra[,152:181])
jul = rowMeans(Ra[,182:212])
aug = rowMeans(Ra[,213:243])
sep = rowMeans(Ra[,244:273])
oct = rowMeans(Ra[,274:304])
nov = rowMeans(Ra[,305:334])
dec = rowMeans(Ra[,335:365])

# Calculate mean daily PE by month
lamda = 2450000
rhow = 1000
Re_Mj = cbind(jan, feb, mar, apr, may, jun, jul, aug, sep, oct, nov, dec)
Re = Re_Mj*1000000
Ta = lat_raw_na[,3:14]

x = Re/(lamda*rhow)
y = (Ta+5)/100
PE_mday = x*y
PE_mmday = PE_mday*1000

# Total PE per month
days = matrix(c(31,28,31,30,31,30,31,31,30,31,30,31), nrow=lat_len, ncol=12, byrow=1)

PE_mmmonth = PE_mmday*days

# Total PE per year
PE_ann = rowSums(PE_mmmonth)

# add x,y columns
PE_monthly = cbind(lat_raw_na[,1:2], PE_mmmonth)
# output results
write.csv(PE_monthly,"~/filepath/PE_monthly_sum.csv",row.names=FALSE)

If you spot a mistake – please let me know!

As a check I plotted up the monthly daily averages for a single location using this code:

# Send result to a png file
png("~/filepath/PEyear.png",width=500,height=400)
# adjust margins
par(mar=c(5,5,2,2)+0.1) # margins: bottom, left, top and right
plot(PE_mmday[1,], xlab="Month of year", ylab=expression(paste("Potential evapotranspiration (mm/day)")), type="l")
dev.off()

Here is the output:

Potential evapotranspiration by month of year for latitude 60.835571.

Potential evapotranspiration by month of year for latitude 60.835571.

I also opened the monthly totals sheet in QGIS and plotted up July:

Monthly total potential evapotranspiration for July.

Monthly total potential evapotranspiration for July.

I’d be really interested to see how this compares to the Met Office MORECS data, so if anyone from the Met Office would like to do a comparison – please get in touch!

The data is available to download here.

About these ads