R Notebook for Water Repellency data process

1 load required libraries and global functions

# check, load and install libraries
if ("dWaSSI" %in% installed.packages()) {library(dWaSSI)} else{
  library(devtools)
  install_github("ln1267/dWaSSI") }
librs<-c("raster","rgdal","rgeos","pryr","RColorBrewer","gstat","rasterVis")
f_lib_check(librs)

2 calculate the index

The band order is: Blue, Green, Red, RedEdge, and Near-Infrared.
\(EG= 2*G-R-B\)
\(NDVI = (NIR - Red)/(NIR+Red)\)
\(NDRE = (NIR - Redge)/(NIR+Redge)\)
\(VARI = (G - R)/(G + R - B)\)

f_cal_VIs<-function(dir_data,Sitename,Tifname,shpname,zonal_field,multiple=T){
  da<-brick(paste0(dir_data,Sitename,"/Original_images/",Tifname))
  shp<-readOGR(paste0(dir_data,Sitename,"/shps"),shpname)
  shp<- spTransform(shp, CRS(proj4string(da)))
  
  f_VARI<-function(red,green,blue){
    
    VARI<-(green-red)/(green+red-blue)
    return(VARI)
  }
  
  f_EG<-function(red,green,blue){
    
    EG<-2*green-red-blue
    return(EG)
  }
  
  f_NDVI<- function(red, nir) {
      ndvi <- (nir - red) / (nir + red)
      return(ndvi)
  }
  f_NDRE<- function(redege, nir) {
      ndre <- (nir - redege) / (nir + redege)
      return(ndre)
  }  
  if(!dir.exists(paste0(dir_data,Sitename,"/Results/Tifs/"))){  dir.create(paste0(dir_data,Sitename,"/Results/Tifs/")) }
  
  print(paste0("Calculating vegetation indices"))

  if(multiple) {
    vis<-c("EG","VARI","NDVI","NDRE")
      # NDVI
    EG <- overlay(da[[3]], da[[2]],da[[1]], fun=f_EG)
    VARI <- overlay(da[[3]], da[[2]],da[[1]], fun=f_VARI)
    NDVI <- overlay(da[[3]],da[[5]], fun=f_NDVI)
    NDRE <- overlay(da[[4]],da[[5]], fun=f_NDRE)
    VARI[is.infinite(VARI) | abs(VARI)>1]<-0
    NDVI[is.infinite(NDVI) | abs(NDVI)>1]<-0
    NDVI[is.infinite(NDRE) | abs(NDRE)>1]<-0
    pdf(paste0(dir_data,Sitename,"/Results/",unlist(strsplit(Tifname, ".", fixed = TRUE))[1],"_NDVI.pdf"),width=14*da@ncols/da@nrows,height=14)
    mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
   plots<-levelplot(NDVI, margin=F, par.settings=mapTheme,main="NDVI",pretty=TRUE)#+layer(sp.points(shp))
   print(plots)
   dev.off()
    #print(plots)
    VIs<-stack(EG,VARI,NDVI,NDRE)
    
  }else{
    vis<-c("EG","VARI")
    # EG
    EG <- overlay(da[[1]], da[[2]],da[[3]], fun=f_EG)
    # vari
    VARI <- overlay(da[[1]], da[[2]],da[[3]], fun=f_VARI)
    VARI[is.infinite(VARI) | abs(VARI)>1]<-0

    pdf(paste0(dir_data,Sitename,"/Results/",unlist(strsplit(Tifname, ".", fixed = TRUE))[1],"_VARI.pdf"),width=14*da@ncols/da@nrows,height=14)
    mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
    plots<-levelplot(VARI, margin=F, par.settings=mapTheme,main="VARI",pretty=TRUE)#+layer(sp.points(shp))
    print(plots)
    dev.off()
    
    VIs<-stack(EG,VARI)
  }

  # zonal data
  print(paste0("Zonaling vegetation index using the field of ",zonal_field))
  ex <- raster::extract(VIs, shp, buffer=0.564,fun=mean,df=T)
  sta_shp <- as.data.frame(ex)
  names(sta_shp)<-c(zonal_field,vis)
  sta_shp[zonal_field] <- shp[[zonal_field]]

  # Export vis
  for( varname in vis){
    print(paste0("exporting ",varname))
    writeRaster(get(varname),paste0(dir_data,Sitename,"/Results/nc/",unlist(strsplit(Tifname, ".", fixed = TRUE))[1],"_",varname,".nc"),varname=varname,xname="lon",yname="lat",format="CDF",overwrite=T)
    writeRaster(get(varname),paste0(dir_data,Sitename,"/Results/Tifs/",unlist(strsplit(Tifname, ".", fixed = TRUE))[1],"_",varname,".tif"),format="GTiff",overwrite=T)
  }
  
  return(sta_shp)
}

f_zonal_VI<-function(dir_data,Sitename,Tifname,shpname,zonal_field,VIname=""){
  da<-raster(paste0(dir_data,Sitename,"/Original_images/",Tifname))
  da[da>1]<-0
  da[da< 0]<-0
  shp<-readOGR(paste0(dir_data,Sitename,"/shps"),shpname)
  shp<- spTransform(shp, CRS(proj4string(da)))
  
  if(!dir.exists(paste0(dir_data,Sitename,"/Results/Tifs/"))){  dir.create(paste0(dir_data,Sitename,"/Results/Tifs/")) }
  
  print(paste0("Zonal vegetation index"))

    pdf(paste0(dir_data,Sitename,"/Results/",unlist(strsplit(Tifname, ".", fixed = TRUE))[1],"_",VIname,".pdf"),width=14*da@ncols/da@nrows,height=14)
    mapTheme <- rasterTheme(region=brewer.pal(8,"RdYlGn"))
    plots<-levelplot(da, margin=F, par.settings=mapTheme,main=VIname,pretty=TRUE)#+layer(sp.points(shp))
    print(plots)
    dev.off()

  # zonal data
  print(paste0("Zonaling vegetation index using the field of ",zonal_field))
  ex <- raster::extract(a, shp, buffer=0.564,fun=mean,df=T)
  sta_shp <- as.data.frame(ex)
  names(sta_shp)<-c(zonal_field,VIname)
  sta_shp[zonal_field] <- shp[[zonal_field]]

  return(sta_shp)
}
#

f_2index_cal<-function(dir_data,filename,daname){
  
  da<-brick(filename)

  # test the input data
  pdf(paste(dir_image,daname,"_RGB.pdf",sep=""))
  #plot(da[[1]])
  plotRGB(da, 1, 2, 3)
  #title(main ="RGB")
  dev.off()
  
  # calculate and write EG
  
  EG <- calc(da, fun=f_EG_raster)
  
  writeRaster(EG,paste(dir_image,daname,"_EG.tif",sep=""),format="GTiff",overwrite=T)
  
  # calculate and write VARI
  VARI<-overlay(da, fun=f_VARI_raster)
  writeRaster(VARI,paste(dir_image,daname,"_VARI.tif",sep=""),format="GTiff",overwrite=T)
  
  # test EG and VARI
  pdf(paste(dir_image,daname,"_index.pdf",sep=""))
  plot(EG)
  title(main ="EG index")
  plot(VARI)
  title(main ="VARI index")
  dev.off()
}


f_2index_array<-function(dir_data,filename,daname,shp=NA){

  da<-brick(filename)
  
  if(!is.na(shp)){
    da<-crop(da,shp)
  }
  
  # test the input data
  pdf(paste(dir_image,daname,"_RGB.pdf",sep=""))
  plotRGB(da, 1, 2, 3,main="RGB image")
  #print(a)
  dev.off()
  
  da_array<-as.array(da)
  # calculate and write EG

  EG <- f_EG_array(da_array)
  EG_raster<-raster(EG,da@extent@xmin,da@extent@xmax,da@extent@ymin,da@extent@ymax,crs=crs(da))
  plot(EG_raster)
  writeRaster(EG_raster,paste(dir_image,daname,"_EG.tif",sep=""),format="GTiff",overwrite=T)

  # calculate and write VARI
  VARI<-f_VARI_array(da_array)
  VARI_raster<-raster(VARI,da@extent@xmin,da@extent@xmax,da@extent@ymin,da@extent@ymax,crs=crs(da))
  writeRaster(VARI_raster,paste(dir_image,daname,"_VARI.tif",sep=""),format="GTiff",overwrite=T)

  # test EG and VARI
  pdf(paste(dir_image,daname,"_index.pdf",sep=""))
  plot(EG_raster)
  title(main ="EG index")
  plot(VARI_raster)
  title(main ="VARI index")
  dev.off()
}

3 Sites

3.1 Meckering

The band order is: Blue, Green, Red, RedEdge, and Near-Infrared.

dir_data<-"/www/Stan/"
# Zonal data for each site
Meckering_VIs_Meckering_trial<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Meckering",
  Tifname="Meckering_trial_site.tif",
  shpname="Meckering_trial",
  zonal_field<-"id",
  multiple=F
  )
save(VIs_Meckering_trial,file = paste0(dir_data,"Tmp/VIs_Meckering_trial.RData"))

Meckering_VIs_082017<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Meckering",
  Tifname="August_multi_2017.tif",
  shpname="Meck_Aug_sample_points",
  zonal_field<-"id"
  )
save(Meckering_VIs_082017,file = paste0(dir_data,"Tmp/Meckering_VIs_082017.RData"))

Meckering_VIs_102017<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Meckering",
  Tifname="October_multi_2017.tif",
  shpname="Meck_Oct_sample_points",
  zonal_field<-"id"
  )
save(Meckering_VIs_102017,file = paste0(dir_data,"Tmp/Meckering_VIs_102017.RData"))

Meckering_VIs_092017<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Meckering",
  Tifname="September_multi_2017.tif",
  shpname="Meck_Sept_sample_points",
  zonal_field<-"id"
  )
save(Meckering_VIs_092017,file = paste0(dir_data,"Tmp/Meckering_VIs_092017.RData"))

3.2 Merge with MED and other data

MED<-read.csv(paste0(dir_data,"Meckering/data/Meckering_med_moist.csv"))
MED<-MED[-10]

load(paste0(dir_data,"Tmp/Meckering_VIs_092017.RData"))
load(paste0(dir_data,"Tmp/Meckering_VIs_082017.RData"))
load(paste0(dir_data,"Tmp/Meckering_VIs_102017.RData"))

Meckering_VIs_092017$Date<-"092017"
Meckering_VIs_082017$Date<-"082017"
Meckering_VIs_102017$Date<-"102017"

VIs<-rbind(Meckering_VIs_092017,Meckering_VIs_082017,Meckering_VIs_102017)
names(VIs)[1]<-"samples"
library(reshape2)
a<-melt(VIs,id=c("samples","Date"))
VIs<-dcast(a,samples~variable+Date)

Merge_data<-merge(MED,VIs,all.x = T,by="samples")

write.csv(Merge_data,paste0(dir_data,Sitename,"/Results/",Sitename,"_Merge_MED_VIs.csv"),row.names = F)

library(ggplot2)
ggplot(data =Merge_data,aes(x=NDVI_092017,y=med_0807,col=Position))+geom_point()+facet_wrap(~Position)

library(dplyr)

cor.test(Merge_data$NDVI,Merge_data$med_0807)

3.3 Kojonup

shpfile<-"Kojonup_2017_sample_points"
# Zonal data for each site
VIs_Kojonup_trial<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Kojonup",
  Tifname="Sept_2017_multi_modified.tif",
  shpname=shpfile,
  zonal_field<-"id",
  multiple=F
  )

shpfile<-"Kojonup_2017_sample_points"
# Zonal data for each site
VIs_Kojonup_trial<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Kojonup",
  Tifname="Kojonup-Tot-Pot-Ura-Tho.tif",
  shpname=shpfile,
  zonal_field<-"id",
  multiple=F
  )

save(VIs_Kojonup_trial,file = "/Dataset/www/Stan/Tmp/VIs_Kojonup_trial.RData")


Kojonup_VIs_092017<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Kojonup",
  Tifname="Sept_2017_multi_modified.tif",
  shpname=shpfile,
  zonal_field<-"id"
  )
save(Kojonup_VIs_092017,file = "/www/Stan/Tmp/Kojonup_VIs_092017.RData")

Kojonup_VIs_102017<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Kojonup",
  Tifname="Oct_2017_multi_modified.tif",
  shpname=shpfile,
  zonal_field<-"id"
  )
save(Kojonup_VIs_102017,file = "/www/Stan/Tmp/Kojonup_VIs_102017.RData")


Kojonup_Flight2<-f_zonal_VI(
  dir_data="/www/Stan/",
  Sitename="Kojonup",
  Tifname="North_Kojonup_Flight_2_NDVI_poly1_lanczos.tif",
  shpname=shpfile,
  zonal_field<-"id",
  VIname = "NDVI"
  )
save(Kojonup_Flight2,file = "/www/Stan/Tmp/Kojonup_Flight2.RData")

plot(a)
plot(shp,add=T)



Kojonup_Flight4<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Kojonup",
  Tifname="North_Kojonup_Flight_4_NDVI_poly1_lanczos.tif",
  shpname=shpfile,
  zonal_field<-"id"
  )
save(Kojonup_Flight4,file = "/www/Stan/Tmp/Kojonup_Flight4.RData")

Kojonup_Flight5<-f_cal_VIs(
  dir_data="/www/Stan/",
  Sitename="Kojonup",
  Tifname="North_Kojonup_Flight_5_NDVI_poly1_lanczos.tif",
  shpname=shpfile,
  zonal_field<-"id"
  )
save(Kojonup_Flight5,file = "/www/Stan/Tmp/Kojonup_Flight5.RData")

3.4 Merge MED

MED<-read.csv("/www/Stan/Kojonup/data/Kojonup_med_moist.csv")

load("/www/Stan/Tmp/Kojonup_VIs_092017.RData")
load("/www/Stan/Tmp/Kojonup_VIs_102017.RData")

Kojonup_VIs_092017$Date<-"092017"
Kojonup_VIs_102017$Date<-"102017"

VIs<-rbind(Kojonup_VIs_092017,Kojonup_VIs_102017)
names(VIs)[1]<-"Samples"
library(reshape2)
a<-melt(VIs,id=c("Samples","Date"))
VIs<-dcast(a,Samples~variable+Date)

Merge_data<-merge(MED,VIs,all.x = T,by="Samples")

write.csv(Merge_data,"/www/Stan/Kojonup/Results/Merge_MED_VIs.csv",row.names = F)

4 Process Radimatrix data

4.1 Sites

# interplate points using IDW 
f_plotrad<-function(da,df,varname){
  library(gstat)
  r<-raster(matrix(1,nrow=as.integer(da@nrows/10),ncol=as.integer(da@ncols/10)),xmn=da@extent@xmin,xmx=da@extent@xmax,ymn=da@extent@ymin,ymx=da@extent@ymax,crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
  # interpolate
  mg <- gstat(id = varname, formula = get(varname)~1, locations = ~x+y, data=df, 
              nmax=7, set=list(idp = .5))
  z <- interpolate(r, mg)
  z <- mask(z, r)
  names(z)<-varname
  mapTheme <- rasterTheme(region=brewer.pal(8,"Blues"))
  plt <- levelplot(z, margin=F, par.settings=mapTheme,main=varname,pretty=TRUE)
  plt
  return(z)
}


da<-raster("/www/Stan/2018/Badgingarra/McAlpinepaddock217072017_VIS.nc")
shp<-readOGR("/www/Stan/2018/Badgingarra","boundary")
print(da)

# read rad EM data and select useful columns
rad<-read.csv("/www/Stan/2018/EM_and_Gamma/EM_and_gamma/EM/Badgingarra EM.csv",header = F,stringsAsFactors=F)


## Load TOTl data
rad<-read.csv("/www/Stan/2018/EM_and_Gamma/EM_and_gamma/Gamma/Badgingarra gamma.csv",header =T,stringsAsFactors=F)
rad<-rad[c(4,3,5:9)]
names(rad)<-c("y","x","Ele","TotCount","Potassium","Uranium","Thorium")

rad<-subset(rad,x>0)
a<-f_plotrad(resolution = 0.00001,df = rad,varname = "TotCount")
b<-f_plotrad(resolution = 0.00001,df = rad,varname = "Potassium")
c<-f_plotrad(resolution = 0.00001,df = rad,varname = "Uranium")
d<-f_plotrad(resolution = 0.00001,df = rad,varname = "Thorium")
cc <- brick(a,b,c,d)

writeRaster(cc,filename = "/www/Stan/2018/EM_and_Gamma/result/Badgingarra-Tot-Pot-Ura-Tho.tif",format="GTiff",overwrite=T)


rad<-read.csv("/www/Stan/2018/EM_and_Gamma/EM_and_gamma/EM/Kojonup EM.csv",header = F,stringsAsFactors=F)

rad<-rad[c(5,6,10,17,18,19,20)]
names(rad)<-c("y","x","Ele","EMd","MG1","EMs","MG2")

rad<-read.csv("/www/Stan/2018/EM_and_Gamma/EM_and_gamma/Gamma/Badgingarra gamma.csv")
rad<-read.csv("/www/Stan/2018/EM_and_Gamma/EM_and_gamma/Gamma/Kojonup gamma.csv")
rad<-rad[3:9]
names(rad)<-c("x","y","Ele","Totalcount","Potassium","Uranium","Thorium")


# transfer csv to spatil points
xy <- rad[,c("x","y")] # long + lat
spdf <- SpatialPointsDataFrame(coords = xy, data = rad,
                               proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

# clip points data based on the outlint polygon
spdf.clip <- spdf[shp, ]

pdf("/www/Stan/gamma.pdf")
for (var in names(rad)[-c(1:3)]){
  
  print(f_plotrad(a,spdf@data,var))

}

dev.off()
for (var in names(rad)[-c(1:3)]){
  
  zz<-f_plotrad(a,spdf@data,var)
  f_plot_sp(zz,filename = paste0("/www/Stan/",var,".pdf"),varnames = var,cuts = 10)
}

r<-f_plotrad(a,spdf@data,"Potassium")
b<-f_plotrad(a,spdf@data,"Uranium")
g<-f_plotrad(a,spdf@data,"Thorium")

c <- brick(r,g,b)
writeRaster(c,"/www/Stan/Kojonup/Results/Kojnup_gamma.tif",format="GTiff",overwrite=T)
pdf("/www/Stan/rad_plot.pdf")
mapTheme <- rasterTheme(region=brewer.pal(8,"Reds"))
levelplot(r, margin=F, par.settings=mapTheme,main="Potassium",pretty=TRUE)
mapTheme <- rasterTheme(region=brewer.pal(8,"Greens"))
levelplot(g, margin=F, par.settings=mapTheme,main="Thorium",pretty=TRUE)
mapTheme <- rasterTheme(region=brewer.pal(8,"Blues"))
levelplot(b, margin=F, par.settings=mapTheme,main="Uranium",pretty=TRUE)

plotRGB(c, 1, 2,3, stretch="hist")

dev.off()
# }

plot(r)
plot(shp,add=T)

4.2 Sites_function

f_gamma(input = "/www/Stan/2018/EM_and_Gamma/Moora_large_scale/DEX_20180212_050745.csv",
       plotname="/www/Stan/2018/EM_and_Gamma/result/DEX-Tot-Pot-Ura-Tho.pdf",
       output="/www/Stan/2018/EM_and_Gamma/result/DEX-Tot-Pot-Ura-Tho.tif",
       plot=T)

[1] "/www/Stan/2018/EM_and_Gamma/Moora_large_scale/DEX_20180212_050745.csv"

Loading required package: lattice
Loading required package: latticeExtra
Loading required package: RColorBrewer

[inverse distance weighted interpolation]
[inverse distance weighted interpolation]
[inverse distance weighted interpolation]
[inverse distance weighted interpolation]