Project

General

Profile

Download (3.6 KB) Statistics
| Branch: | Revision:
1
# R code to plot latitudinal profiles of mean slope, comparing fused layers 
2
# with both the raw ASTER. I did not include the comparison of both RMSE 
3
# and correlation coefficients,which Jim included in his analysis. 
4
#
5
# The results (plots) are saved at 
6
#   "/data/project/organisms/DEM?Yuni/documents/slope"
7
#
8
# Original author: Jim Regetz
9
# [13-JUL-2011]
10
#   Edits by Yuina to auto-apply for nine regions.
11
#   
12
#   9-Dec-2011
13
#   Yuina Nunokawa 
14

    
15

    
16

    
17

    
18
library(raster)
19

    
20

    
21
# load slope rasters
22
aster.dir      <- "/data/project/organisms/DEM/Yuni/Data/aster2"
23
srtm.dir       <- "/data/project/organisms/DEM/Yuni/Data/srtm"
24
aster.straddle <- list.files(aster.dir, pattern="^aster2.*_straddle_s.tif")
25
srtm.below     <- list.files(srtm.dir, pattern="^srtm_.*_below_s.tif") 
26
fused.straddle <- list.files(aster.dir, pattern="^fused.*_straddle_s.tif")
27
fused.bgau     <- list.files(aster.dir, pattern="^fused.*_blendgau_s.tif")
28

    
29

    
30
slope <- function(s.aster2, s.srtm, s.uncor, s.bg, file.name, ylim, aster.dir, srtm.dir){
31

    
32
        # extract raster latitudes for later
33
        lats2400 <- yFromRow(s.aster2, 1:nrow(s.aster2))
34
        lats1200 <- yFromRow(s.srtm, 1:nrow(s.srtm))
35

    
36
        # initialize output pdf device driver
37
        name <- paste("slope_", file.name, ".pdf", sep="") 
38
        pdf(file.path("/data/project/organisms/DEM/Yuni/documents/slope", file=name), height=8, width=11.5)
39

    
40
        #
41
        # plot latitudinal profiles of mean slope
42
        #
43
        par(mfrow=c(2,2), omi=c(1,1,1,1))
44
	ylim <- ylim
45
	plot(lats2400, rowMeans(as.matrix(s.aster2), na.rm=TRUE), type="l", xlab="Latitude", ylab="Mean slope", ylim=ylim)
46
        lines(lats1200, rowMeans(as.matrix(s.srtm), na.rm=TRUE), col="red")
47
        text(min(lats2400), max(ylim)-0.5, pos=4, font=3, labels="Original")
48
        legend("bottomright", legend=c("ASTER2", "SRTM"), col=c("blue", "red"), lty=c(1, 1), bty="n")
49
        abline(v=60, col="red", lty=2)
50
        mtext(expression(paste("Latitudinal profiles of mean slope", file.name)), adj=0, line=2, font=2)
51

    
52
        plot(lats2400, rowMeans(as.matrix(s.uncor), na.rm=TRUE), type="l", xlab="Latitude", ylab="Mean slope", ylim=ylim)
53
        legend("bottomright", legend=c("ASTER2"), col=c("black"), lty=c(1, 1), bty="n")
54
        text(min(lats2400), max(ylim)-0.5, pos=4, font=3, labels="simple fused")
55
        abline(v=60, col="red", lty=2)
56

    
57
        plot(lats2400, rowMeans(as.matrix(s.bg), na.rm=TRUE), type="l", xlab="Latitude", ylab="Mean slope", ylim=ylim)
58
        legend("bottomright", legend=c("ASTER2"), col=c("black"), lty=c(1, 1), bty="n")
59
        text(min(lats2400), max(ylim)-0.5, pos=4, font=3, labels="gaussian blend")
60
        abline(v=60, col="red", lty=2)
61
    
62
        dev.off () 
63

    
64
}
65

    
66

    
67

    
68

    
69
num.regions <- 9
70
plotSlope <-lapply(1:num.regions, function(region, aster.dir, srtm.dir){
71
    s.aster2 <- raster(file.path(aster.dir, aster.straddle[region])) 
72
    s.srtm <- raster(file.path(srtm.dir, srtm.below[region]))  
73
    s.uncor <- raster(file.path(aster.dir, fused.straddle[region]))
74
    s.bg <- raster(file.path(aster.dir, fused.bgau[region])) 
75
    file.name <- substr(srtm.below[region], 6, 13)
76
    
77
    if (region==1){
78
        ylim <- c(0, 4)
79
    }else if(region==2){
80
        ylim <- c(0, 4)
81
    }else if (region==3){   
82
        ylim <- c(2, 7)
83
    }else if (region==4){
84
        ylim <-c(1, 7)
85
    }else if (region==5){
86
        ylim <- c(0, 3)
87
    }else if (region==6){
88
        ylim <- c(0, 8)
89
    }else if (region==7){
90
        ylim <- c(0, 3)
91
    }else if (region==8){
92
        ylim <- c(3, 7)
93
    }else{
94
        ylim <- c(0, 6)
95
    }
96
    slope(s.aster2, s.srtm, s.uncor, s.bg, file.name, ylim, aster.dir, srtm.dir)
97
}, aster.dir=aster.dir, srtm.dir=srtm.dir)
98

    
99

    
(22-22/23)