I am currently helping supervise a masters student for their summer project. They’re working on improving a parameter used in flood estimation. The drive for the project is to use the revised parameter in the assessment of natural flood management, something SEPA and local authorities are now bound to consider since the implementation of the Flood Risk Management (Scotland) Act 2009.
The parameter derivation is dependant on a number of datasets and uses these to undertake a spatial regression, thus offering a new definition.
One dataset is a grid of potential evapotranspiration. For national assessments in the UK this is traditionally taken from the MORECS data owned by the Met Office. Unfortunately the cost of this was prohibitive (£4000 plus VAT for the UK) and as a core product, the Met Office were unwilling to release it with a student licence.
With a bit of hunting I came across this article by Kay and Davies, which derives a PE grid from freely available data and makes comparisons to MORECS. This seemed a viable solution and so became our chosen path.
To calculate the PE data extraterrestrial radiation is required. Kay and Davies’ paper provides a reference to Allen et al. who calculates this based on Julian day and latitude (calcs on pdf page 24, journal page 70). As this has potential use for my work (input solar energy probably has some impact on snow accumulation and melt!) and the student was struggling: I decided to roll up my sleeves and do some sums. You can see page 70 from the paper below, showing the required calculations.
- Created a 5 km grid to match the Met Office UKCIP09 one, because we’ll be using the 5 km temperature grid to calculate potential evapotranspiration (using QGIS Vector/Research Tools/Vector grid)
- Extracted the polygon centroids from the vector grid (QGIS Vector/Geometry Tools/Polygon centroids)
- Previous two steps were working in British National Grid, EPSG:27700. I converted the layer to WGS84 (EPSG:4326), as I needed latitude, not northing (QGIS save as…)
- Opened the new WGS84 file in QGIS and appended geometry columns (QGIS Vector/Geometry Tools/Export/Add geometry columns)
- Exported the file as a csv (QGIS save as…)
I next spent a bit of time writing a script in R to do the work. Initially I toyed with the idea of using a spreadsheet, but the thought of many fill downs and all the potential errors were not appealing. Most of the script writing time was spent working out how to get the matrices to appear in the sums correctly.
If you’re reading this and you think I’ve made a mistake – please tell me!
# import latitude lat_raw = as.matrix(read.csv("~/file path/ext_rad_wgs84_cent.csv",sep=",",header=1)) # prepare latitude lat_len = nrow(lat_raw) lat_rad = lat_raw[ ,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 extraterrestrial radiation 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)) # sum every day to get annual Ra_ann = rowSums(Ra) # add x,y columns Ra_ann_output = cbind(lat_raw[,1:2], Ra_ann) # output results write.csv(Ra_ann_output,"~/file path/rad_ann.csv",row.names=0)
The script outputs to a csv file with two columns of coordinates (lat-long in degrees) with a third column of annual extraterrestrial radiation. As a quick check I plotted up the annual cycle of radiation for one latitude, expecting a bell curve (or sine curve, depending on how you look at it). Which I got.
I also plotted the annual data on British National Grid, EPSG:27700 and colour coded by total radiation. As expected, further north gets less radiation.
If you’d like the data, please download it here.