Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # calculate the warmth index
- calcWI = function(dir, gcm, rcp, epoch, filetype='.tif') {
- tx = list(); tn = list()
- tm = list(); tm_adj = list()
- for ( i in 1:12 ) {
- tx[[i]] = raster::raster(rgdal::readGDAL(paste(dir, 'tw_', gcm, rcp,
- 'tx', epoch, i, filetype, sep='')))/10.0
- tn[[i]] = raster::raster(rgdal::readGDAL(paste(dir, 'tw_', gcm, rcp,
- 'tn', epoch, i, filetype, sep='')))/10.0
- tm[[i]] = ( tx[[i]] + tn[[i]] ) / 2.0
- tm_adj[[i]] = tm[[i]]-5
- tm_adj[[i]][tm_adj[[i]]<0] = 0
- if ( i == 1 ) {
- wi = tm_adj[[i]]
- } else {
- wi = wi + tm_adj[[i]]
- }
- opt = paste(dir, 'tw_', gcm, rcp, 'wi', epoch, '.asc', sep='')
- raster::writeRaster(wi, opt, format='ascii',overwrite=TRUE)
- }
- }
- gcm=c('bc', 'cc', 'gs', 'hd', 'he', 'ip', 'mc', 'mg', 'mi', 'mr', 'no')
- rcp=c('26','45','60','85')
- epoch=c('50','70')
- pre='/Volumes/Data/cmip5'
- for ( i in 1:length(gcm) ) {
- for ( j in 1:length(rcp) ) {
- for ( k in 1:length(epoch) ) {
- gcm_dir = paste(pre,gcm[i],'/',sep='/')
- calcWI(dir=gcm_dir, gcm=gcm[i], rcp=rcp[j], epoch=epoch[k], filetype='.tif')
- }
- }
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement