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.
- 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).
# 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:
I also opened the monthly totals sheet in QGIS and plotted up 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.