Advertisement
Guest User

Untitled

a guest
Apr 1st, 2015
222
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.18 KB | None | 0 0
  1. # calculate the warmth index
  2. calcWI = function(dir, gcm, rcp, epoch, filetype='.tif') {
  3. tx = list(); tn = list()
  4. tm = list(); tm_adj = list()
  5. for ( i in 1:12 ) {
  6. tx[[i]] = raster::raster(rgdal::readGDAL(paste(dir, 'tw_', gcm, rcp,
  7. 'tx', epoch, i, filetype, sep='')))/10.0
  8. tn[[i]] = raster::raster(rgdal::readGDAL(paste(dir, 'tw_', gcm, rcp,
  9. 'tn', epoch, i, filetype, sep='')))/10.0
  10. tm[[i]] = ( tx[[i]] + tn[[i]] ) / 2.0
  11. tm_adj[[i]] = tm[[i]]-5
  12. tm_adj[[i]][tm_adj[[i]]<0] = 0
  13. if ( i == 1 ) {
  14. wi = tm_adj[[i]]
  15. } else {
  16. wi = wi + tm_adj[[i]]
  17. }
  18. opt = paste(dir, 'tw_', gcm, rcp, 'wi', epoch, '.asc', sep='')
  19. raster::writeRaster(wi, opt, format='ascii',overwrite=TRUE)
  20. }
  21. }
  22. gcm=c('bc', 'cc', 'gs', 'hd', 'he', 'ip', 'mc', 'mg', 'mi', 'mr', 'no')
  23. rcp=c('26','45','60','85')
  24. epoch=c('50','70')
  25. pre='/Volumes/Data/cmip5'
  26. for ( i in 1:length(gcm) ) {
  27. for ( j in 1:length(rcp) ) {
  28. for ( k in 1:length(epoch) ) {
  29. gcm_dir = paste(pre,gcm[i],'/',sep='/')
  30. calcWI(dir=gcm_dir, gcm=gcm[i], rcp=rcp[j], epoch=epoch[k], filetype='.tif')
  31. }
  32. }
  33. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement