Load data

load("data/GewataB1.rda")
load("data/GewataB2.rda")
load("data/GewataB3.rda")
load("data/GewataB4.rda")
load("data/GewataB5.rda")
load("data/GewataB7.rda")
load("data/vcfGewata.rda")

Clean water, clouds and cloud shadow from VCF data

vcfGewata[vcfGewata>100]<-NA
hist(vcfGewata)

Build a brick containing all data

gwt <- brick(GewataB1, GewataB2, GewataB3, GewataB4, GewataB5, GewataB7, vcfGewata)
names(gwt) <- c("band1", "band2", "band3", "band4", "band5", "band7", "vcf")

Extract all data to a data.frame

gwt.df <- as.data.frame(getValues(gwt)) 
head(gwt.df)

##   band1 band2 band3 band4 band5 band7 vcf
## 1   376   519   520  2479  1819  1102  49
## 2   397   563   577  2310  1868  1153  37
## 3   439   585   674  2140  1966  1254  28
## 4   439   563   694  2054  1600  1026  28
## 5   439   651   771  2096  1844  1331  26
## 6   501   739   943  2351  2379  1610  19

str(gwt.df)

## 'data.frame':    1821996 obs. of  7 variables:
##  $ band1: int  376 397 439 439 439 501 521 501 439 419 ...
##  $ band2: int  519 563 585 563 651 739 739 673 629 563 ...
##  $ band3: int  520 577 674 694 771 943 905 789 712 615 ...
##  $ band4: int  2479 2310 2140 2054 2096 2351 2310 2181 2181 2225 ...
##  $ band5: int  1819 1868 1966 1600 1844 2379 2403 2209 1990 1795 ...
##  $ band7: int  1102 1153 1254 1026 1331 1610 1712 1534 1331 1153 ...
##  $ vcf  : int  49 37 28 28 26 19 16 18 28 25 ...

Assign bands and VCF values to vectors

band1<-gwt.df$band1
band2<-gwt.df$band2
band3<-gwt.df$band3
band4<-gwt.df$band4
band5<-gwt.df$band5
band7<-gwt.df$band7
vcf<-gwt.df$vcf

Visualizing the data and the relationship between bands and VCF tree cover

plot(gwt)

summary(gwt)

##         gewataB1 gewataB2 gewataB3 gewataB4 gewataB5 gewataB7
## Min.          89      127       85      598      113       34
## 1st Qu.      326      410      307     2178     1139      491
## Median       373      469      368     2407     1357      618
## 3rd Qu.      441      594      582     2646     1808     1051
## Max.        1545     2220     2865     4903    16000    16000
## NA's           0        0        0        0        0        0
##         vcf2000Gewata
## Min.                0
## 1st Qu.            32
## Median             64
## 3rd Qu.            75
## Max.              100
## NA's            13712

hist(gwt)
pairs(gwt)

Creat a lm model

lm1<-lm(vcf~band1+band2+band3+band4+band5+band7)
summary(lm1)

## 
## Call:
## lm(formula = vcf ~ band1 + band2 + band3 + band4 + band5 + band7)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -58.205  -4.636   0.715   5.211 258.436 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  8.549e+01  5.723e-02 1493.805   <2e-16 ***
## band1        9.291e-02  1.889e-04  491.817   <2e-16 ***
## band2       -1.505e-01  2.353e-04 -639.622   <2e-16 ***
## band3       -2.528e-03  1.698e-04  -14.889   <2e-16 ***
## band4        1.661e-02  2.790e-05  595.378   <2e-16 ***
## band5       -2.006e-02  7.036e-05 -285.119   <2e-16 ***
## band7        2.549e-05  8.969e-05    0.284    0.776    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.606 on 1808277 degrees of freedom
##   (13712 observations deleted due to missingness)
## Multiple R-squared:  0.8582, Adjusted R-squared:  0.8582 
## F-statistic: 1.824e+06 on 6 and 1808277 DF,  p-value: < 2.2e-16

From the results, we can see that band 7 is of little importance

in predicting VCF. For instance, if we remove band 7, the resulting

R-squared would still be the same:

summary(lm(vcf~band1+band2+band3+band4+band5))

## 
## Call:
## lm(formula = vcf ~ band1 + band2 + band3 + band4 + band5)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -58.204  -4.636   0.715   5.211 258.586 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.550e+01  5.601e-02  1526.5   <2e-16 ***
## band1        9.291e-02  1.876e-04   495.3   <2e-16 ***
## band2       -1.505e-01  2.314e-04  -650.5   <2e-16 ***
## band3       -2.504e-03  1.473e-04   -17.0   <2e-16 ***
## band4        1.661e-02  2.443e-05   679.9   <2e-16 ***
## band5       -2.004e-02  4.158e-05  -482.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.606 on 1808278 degrees of freedom
##   (13712 observations deleted due to missingness)
## Multiple R-squared:  0.8582, Adjusted R-squared:  0.8582 
## F-statistic: 2.189e+06 on 5 and 1808278 DF,  p-value: < 2.2e-16

Predict tree cover from model

predTC<- predict(gwt,model=lm1,na.rm=TRUE)
par(mfrow=c(1,1))
predTC[predTC<0]<-NA
plot(predTC,xlim=c(808755,855195),ylim=c(817635,852945))

str(predTC)

## Formal class 'RasterLayer' [package "raster"] with 12 slots
##   ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. ..@ name        : chr ""
##   .. .. ..@ datanotation: chr "FLT4S"
##   .. .. ..@ byteorder   : chr "little"
##   .. .. ..@ nodatavalue : num -Inf
##   .. .. ..@ NAchanged   : logi FALSE
##   .. .. ..@ nbands      : int 1
##   .. .. ..@ bandorder   : chr "BIL"
##   .. .. ..@ offset      : int 0
##   .. .. ..@ toptobottom : logi TRUE
##   .. .. ..@ blockrows   : int 0
##   .. .. ..@ blockcols   : int 0
##   .. .. ..@ driver      : chr ""
##   .. .. ..@ open        : logi FALSE
##   ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. ..@ values    : num [1:1821996] 45.7 37.1 32.7 41.8 24.2 ...
##   .. .. ..@ offset    : num 0
##   .. .. ..@ gain      : num 1
##   .. .. ..@ inmemory  : logi TRUE
##   .. .. ..@ fromdisk  : logi FALSE
##   .. .. ..@ isfactor  : logi FALSE
##   .. .. ..@ attributes: list()
##   .. .. ..@ haveminmax: logi TRUE
##   .. .. ..@ min       : num 0.000956
##   .. .. ..@ max       : num 97.9
##   .. .. ..@ band      : int 1
##   .. .. ..@ unit      : chr ""
##   .. .. ..@ names     : chr ""
##   ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. ..@ type      : chr(0) 
##   .. .. ..@ values    : logi(0) 
##   .. .. ..@ color     : logi(0) 
##   .. .. ..@ names     : logi(0) 
##   .. .. ..@ colortable: logi(0) 
##   ..@ title   : chr(0) 
##   ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. ..@ xmin: num 808755
##   .. .. ..@ xmax: num 855195
##   .. .. ..@ ymin: num 817635
##   .. .. ..@ ymax: num 852945
##   ..@ rotated : logi FALSE
##   ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. ..@ geotrans: num(0) 
##   .. .. ..@ transfun:function ()  
##   ..@ ncols   : int 1548
##   ..@ nrows   : int 1177
##   ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=utm +zone=36 +ellps=WGS84 +units=m +no_defs"
##   ..@ history : list()
##   ..@ z       : list()

randomForest model

## Load the calibration points
samples<-readOGR("data/calib.gwt.kml",layer="calib.gwt")

## OGR data source with driver: KML 
## Source: "data/calib.gwt.kml", layer: "calib.gwt"
## with 16 features
## It has 2 fields

class(samples)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

## Re-project spatialPointsDataFrame
samples2<-spTransform(samples,CRS(proj4string(gwt)))
samples2@coords<-coordinates(samples2)[,-3]
## Extract the surface reflectance
calib<- extract(gwt,samples2,df=TRUE)
## Combine the created DF to the description columm of the calibration dataset
calib2<- cbind(samples2$Description,calib)
calib2<-na.omit(calib2)
str(calib2)

## 'data.frame':    13 obs. of  9 variables:
##  $ samples2$Description: Factor w/ 3 levels "crop","forest",..: 2 2 1 1 1 2 2 3 3 2 ...
##  $ ID                  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ band1               : int  475 380 565 530 593 424 396 458 546 382 ...
##  $ band2               : int  517 489 786 677 737 489 501 613 710 414 ...
##  $ band3               : int  393 307 784 748 874 347 378 516 696 339 ...
##  $ band4               : int  3027 2685 2364 2274 1923 2946 2632 2764 2352 2811 ...
##  $ band5               : int  1414 1267 2709 2486 2160 1416 1240 1537 1943 1453 ...
##  $ band7               : int  595 518 1773 1565 1610 570 595 724 1052 592 ...
##  $ vcf                 : int  70 65 12 14 25 70 67 33 29 68 ...
##  - attr(*, "na.action")=Class 'omit'  Named int [1:3] 13 14 15
##   .. ..- attr(*, "names")= chr [1:3] "13" "14" "15"

Calibrate a random forest model using the extracted data frame

rFm<- randomForest(vcf~band1+band2+band3+band4+band5+band7,data=calib2)

Use the model to predict tree cover

tcMap<-predict(gwt,model=rFm)
plot(tcMap)

Function to compute the RMSE

RMSE<-function(x,y) {
  sqrt(mean((y-x$vcf)^2,na.rm=TRUE))
}

## Compute the RMSE for the lm model
predTC.df<-as.data.frame(getValues(predTC))
RMSE(gwt.df,predTC.df)

## [1] 8.356825

## Compute the RMSE for the random forest model
tcMap.df<-as.data.frame(getValues(tcMap))
RMSE(gwt.df,tcMap.df)

## [1] 9.259979

Compute RMSE for each class

## Make an NA-value raster based on the LC raster attributes
formask <- setValues(raster(samples2), NA)
## Assign 1 to formask to all cells corresponding to the forest class
formask[samples2$Description=="forest"] <- 1
## Create a grass mask
grassmask <- setValues(raster(samples2), NA)
grassmask[samples2$Description=="grass"] <- 1
##Create a crop mask
cropmask <- setValues(raster(samples2), NA)
cropmask[samples2$Description=="crop"] <- 1

Get values for each class

tcfor<-crop(tcMap,formask)
tcfor.df<-as.data.frame(getValues(tcfor))
tcgrass<-crop(tcMap,grassmask)
tcgrass.df<-as.data.frame(getValues(tcgrass))
tccrop<-crop(tcMap,cropmask)
tccrop.df<-as.data.frame(getValues(tccrop))

Compute RMSE for each class

RMSE(gwt.df,tcfor.df)

## [1] 31.74691

RMSE(gwt.df,tcgrass.df)

## [1] 31.74691

RMSE(gwt.df,tccrop.df)

## [1] 31.74691

Exercise 8

Paulo Bernardino

4 de fevereiro de 2017

Load libraries

Load data

Clean water, clouds and cloud shadow from VCF data

Build a brick containing all data

Extract all data to a data.frame

Assign bands and VCF values to vectors

Visualizing the data and the relationship between bands and VCF tree cover

Creat a lm model

From the results, we can see that band 7 is of little importance

in predicting VCF. For instance, if we remove band 7, the resulting

R-squared would still be the same:

Predict tree cover from model

randomForest model

Calibrate a random forest model using the extracted data frame

Use the model to predict tree cover

Function to compute the RMSE

Compute RMSE for each class

Get values for each class

Compute RMSE for each class