The relationship ndvi vs. savanna tree density: the segmentation approach

knitr::opts_chunk$set(fig.path = "figure/RmcgwireNewv20180417_")

Preliminaries

#rwd <- "/media/Iomega_HDD/mcgwire/Rmcgwire"
rastdir <- ".."
vectdir <- "../OTB"

require(rgdal)
require(raster)
#require(gplots)
require(ggplot2)
require(gridExtra)
require(colorRamps)
require(RStoolbox)
source("mcgwirecoarsen.R")

(continuation of RmcgwireNewv20180417a_log.html) ## 5. The segmentation alternative

segm <- readOGR(dsn=vectdir, layer="tm_all_step3",stringsAsFactors=FALSE)
## OGR data source with driver: ESRI Shapefile 
## Source: "../OTB", layer: "tm_all_step3"
## with 186 features
## It has 16 fields
## Integer64 fields read as strings:  label nbPixels
#NOTE: try improving these plots with ggplot, see https://github.com/tidyverse/ggplot2/wiki/plotting-polygon-shapefiles
par(mfrow=c(1,2))
plotRGB(tm,r=2,g=4,b=3,stretch="lin")
plot(segm,add=TRUE)
title(main="segm10")
segm10 <- readOGR(dsn=vectdir, layer="tm_all_step3_10",stringsAsFactors=FALSE)
## OGR data source with driver: ESRI Shapefile 
## Source: "../OTB", layer: "tm_all_step3_10"
## with 744 features
## It has 16 fields
## Integer64 fields read as strings:  label nbPixels
plotRGB(tm,r=2,g=4,b=3,stretch="lin")
plot(segm10,add=TRUE)
title(main="segm10")

5.1 Calculate ndv vs. Tree density relationship by segments

#avoid next 2 lines in interactive (slow)
#segm10trees <- extract(trees10,segm10,df=TRUE,fun=mean,na.rm=TRUE)
#segm10ndvi <-  extract(ndvi,segm10, df=TRUE,fun=mean,na.rm=TRUE)
#segmtrees <- extract(trees10,segm,df=TRUE,fun=mean,na.rm=TRUE)
#segmndvi <-  extract(ndvi,segm, df=TRUE,fun=mean,na.rm=TRUE)

segm10ndvitrees <- data.frame(segm10ndvi,percenttrees=segm10trees[,2])
segm10ndvitrees[1:3,]
##   ID      ndvi percenttrees
## 1  1 0.1957324     42.85185
## 2  2 0.1258014     13.53846
## 3  3 0.1387881     14.66667
ggplot(data=segm10ndvitrees,aes(x=percenttrees, y=ndvi)) +
  geom_bin2d(bins=50) +
  scale_fill_gradientn(colours=c("grey70",matlab.like(9))) +
  geom_smooth(method='lm', col="black") +
  ggtitle("Finer Segmentation")