#This script processes the MODIS data with Hansen tree cover

library(rgdal)
library(tiff)
library(raster)
library(ncdf4)

threshold<-80 #minimun mean forest cover in Modis pixel

site.names<-c("Hainich","Bialowieza","Dundo","Kalkalpen","Laegeren","PenedaGeres","Seitseminen")
filenames<-c("Hansen_GFC-2016-v1.4_treecover2000_60N_010E.tif",
              "Hansen_GFC-2016-v1.4_treecover2000_60N_020E.tif",
              "Hansen_GFC-2016-v1.4_treecover2000_50N_010E.tif",
              "Hansen_GFC-2016-v1.4_treecover2000_50N_010E.tif",
              "Hansen_GFC-2016-v1.4_treecover2000_50N_000E.tif",
              "Hansen_GFC-2016-v1.4_treecover2000_50N_010W.tif",
              "Hansen_GFC-2016-v1.4_treecover2000_70N_020E.tif")

path.hansen<-"D:/Maik/Documents/Data/Hansen_forest_cover/"
ls<-list()
for(i in 1:length(site.names)){
  ### load Modis data
  path<-paste("D:/Maik/Documents/Data/MODIS/data/nc/",site.names[i],"/MOD17A3H.006_aid0001.nc",sep="")
  r_modis <- raster(path,  varname = "Npp_500m")
  proj4string(r_modis)=CRS("+init=EPSG:4326")
  
  plot(r_modis)
  
  ### load Hansen data
  Rfd<-readGDAL(paste(path.hansen,filenames[i],sep="")) # read tiff
  r_hansen_large<-raster(Rfd)
  
  r_hansen<-r_modis # copy modis raster
  values(r_hansen)<-1 # create modis mask (every cell values is 1)
  
  hansen_cells<-extract(r_hansen_large,polygons(rasterToPolygons(r_hansen,dissolve = F))) #extract values of hansen raster, masked by modis raster cells
  mean_hansen_cells<-lapply(hansen_cells,FUN=mean,na.rm=T) #aggregate hansen cell values
  
  for(icell in 1:length(r_hansen)){r_hansen[icell]<-mean_hansen_cells[[icell]]} #write aggregated values to raster
  
  r_hansen[is.na(r_modis)]<-NA #replace originally empty Modis pixel with NA
  
  #mask for cutting modis raster:
  r_modis_cut<-r_hansen>threshold
  #if(i==7){ r_modis_cut<-r_hansen>75}
  
  ls[[i]]<-Which(r_modis_cut == 1, cells = TRUE) 
  
  
  cut_hansen<-crop(r_hansen_large,r_modis)
  plot(cut_hansen)
  plot(r_hansen)
  plot(r_modis_cut)
  
}


for(i in 1:length(site.names)){
  write.csv(ls[[i]],file=paste("D:/Maik/Documents/Data/MODIS/data/nc/",site.names[i],"_Modispixels_80.csv",sep=""),row.names = F,col.names = F)
}