############################################################################ # # Author: Daniel M. Griffith, Stephanie Pau, and Ryan Slapikas # Date: 1-June-2020 modified 9-Feb-2021 # Description: Download AOP data, clip to Konza, create NDVI mask # ############################################################################ # SET WD TO SCRIPT LOCATION #setwd(".") # LOAD LIBRARIES library(devtools) library(neonUtilities) library(raster) #devtools::install_github("NEONScience/NEON-geolocation/geoNEON") library(geoNEON) #install.packages("BiocManager") #BiocManager::install("rhdf5") library(rhdf5) library(rgdal) library(rgl) library(rasterVis) library(plyr) #library(maptools) library(foreach) library(doSNOW) # Set global option to NOT convert all character variables to factors options(stringsAsFactors=F) dir.create(path = "neon-downloads", showWarnings = FALSE) # kz site # Meta sitebnds <- readOGR(dsn = "fieldSamplingBoundaries", layer = "terrestrialSamplingBoundaries") # kz # requires you to enter y/n so running this as a full script skips over this part zipsByProduct(dpID = "DP1.10058.001", site = c("KONZ"), package = "expanded", startdate = NA, enddate = NA, # indicates all dates available savepath = "neon-downloads") # neon-downloads/filesToStack10058 stackByTable(filepath = "neon-downloads/filesToStack10058", savepath = "neon-out", saveUnzippedFiles = T, nCores = 1) # Lidar # DP3.30015.001 dir.create(path = "neon-lidar", showWarnings = FALSE) byFileAOP(dpID = "DP3.30015.001", site = "KONZ", year = "2017", savepath = "neon-lidar", check.size = F) # AOP nit - this product has been suspended by NEON and data are only available upon request # https://urldefense.com/v3/__https://data.neonscience.org/data-products/DP3.30018.001__;!!PhOWcWs!mmQ-NjPfqPvsswEhj8X-JePbfMpp68oEteMHJ5NcMOZ759Oq7vP3QKR70lOF$ # DP3.30018.001 dir.create(path = "neon-nit", showWarnings = FALSE) byFileAOP(dpID = "DP3.30018.001", site = "KONZ", year = "2017", savepath = "neon-nit", check.size = F) # Elevation # DP3.30024.001 dir.create(path = "neon-elv", showWarnings = FALSE) byFileAOP(dpID = "DP3.30024.001", site = "KONZ", year = "2017", savepath = "neon-elv", check.size = F) # LAI # DP3.30012.001 dir.create(path = "neon-lai", showWarnings = FALSE) byFileAOP(dpID = "DP3.30012.001", site = "KONZ", year = "2017", savepath = "neon-lai", check.size = F) # BIO - this product has been suspended by NEON and data are only available upon request # https://urldefense.com/v3/__https://data.neonscience.org/data-products/DP3.30016.001__;!!PhOWcWs!mmQ-NjPfqPvsswEhj8X-JePbfMpp68oEteMHJ5NcMOZ759Oq7vP3QAyt9Rbm$ # DP3.30016.001 dir.create(path = "neon-bio", showWarnings = FALSE) byFileAOP(dpID = "DP3.30016.001", site = "KONZ", year = "2017", savepath = "neon-bio", check.size = F) # full spec # DP3.30006.001 dir.create(path = "neon-spc", showWarnings = FALSE) byFileAOP(dpID = "DP3.30006.001", site = "KONZ", year = "2017", savepath = "neon-spc", check.size = F) byTileAOP(dpID="DP3.30006.001", site="KONZ", year="2017", easting=c(717000,717000), northing=c(4335000,4335000)) ################################################################################ ################################################################################ ################################################################################ # multiple flight paths; mosaicked; chm = canopy height model; each file is a flight path; function raster mosaics them; any overlap will choose the max value dir.create(path = "KONZ_mosaics", showWarnings = FALSE) sitebnds_KONZ <- sitebnds[sitebnds$siteID == "KONZ",] chm_files <- list.files(path = "neon-lidar/DP3.30015.001/2017/FullSite/D06/2017_KONZ_2/L3/DiscreteLidar/CanopyHeightModelGtif/", full.names = T, pattern = ".tif") chm_rasters <- lapply(X = chm_files, FUN = raster) chm_rasters$fun <- max chm_rasters$na.rm <- TRUE chm <- do.call(what = mosaic, args = chm_rasters) # used later for other products sitebnds_KONZ <- spTransform(x = sitebnds_KONZ, CRSobj = crs(chm)) chm <- mask(x = chm, mask = sitebnds_KONZ) chm <- crop(x = chm, y = sitebnds_KONZ) writeRaster(x = chm, filename = "KONZ_mosaics/KONZ_chm.tif") ###### dsm dsm_files <- list.files(path = "neon-elv/DP3.30024.001/2017/FullSite/D06/2017_KONZ_2/L3/DiscreteLidar/DSMGtif/", full.names = T, pattern = ".tif") dsm_rasters <- lapply(X = dsm_files, FUN = raster) dsm_rasters$fun <- max dsm_rasters$na.rm <- TRUE dsm <- do.call(what = mosaic, args = dsm_rasters) dsm <- mask(x = dsm, mask = sitebnds_KONZ) dsm <- crop(x = dsm, y = sitebnds_KONZ) writeRaster(x = dsm, filename = "KONZ_mosaics/KONZ_dsm.tif") ###### ndni # Veg Indices in zipped files nit_files <- list.files(path = "neon-nit/DP3.30018.001/2017/FullSite/D06/2017_KONZ_2/L3/Spectrometer/VegIndices/", full.names = T, recursive = T, pattern = "NDNI.tif") nit_rasters <- lapply(X = nit_files, FUN = raster) nit_rasters$fun <- max nit_rasters$na.rm <- TRUE nit <- do.call(what = mosaic, args = nit_rasters) nit <- mask(x = nit, mask = sitebnds_KONZ) nit <- crop(x = nit, y = sitebnds_KONZ) writeRaster(x = nit, filename = "KONZ_mosaics/KONZ_nit.tif") ###### lai lai_files <- list.files(path = "neon-lai/DP3.30012.001/2017/FullSite/D06/2017_KONZ_2/L3/Spectrometer/LAI/", full.names = T, recursive = T, pattern = "LAI.tif") lai_rasters <- lapply(X = lai_files, FUN = raster) lai_rasters$fun <- max lai_rasters$na.rm <- TRUE lai <- do.call(what = mosaic, args = lai_rasters) lai <- mask(x = lai, mask = sitebnds_KONZ) lai <- crop(x = lai, y = sitebnds_KONZ) writeRaster(x = lai, filename = "KONZ_mosaics/KONZ_lai.tif") ###### bio bio_files <- list.files(path = "neon-bio/DP3.30016.001/2017/FullSite/D06/2017_KONZ_2/L3/Spectrometer/Biomass/", full.names = T, recursive = T, pattern = "mass.tif") bio_rasters <- lapply(X = bio_files, FUN = raster) bio_rasters$fun <- max bio_rasters$na.rm <- TRUE bio <- do.call(what = mosaic, args = bio_rasters) bio <- mask(x = bio, mask = sitebnds_KONZ) bio <- crop(x = bio, y = sitebnds_KONZ) writeRaster(x = bio, filename = "KONZ_mosaics/KONZ_bio.tif") # NDVI mask ndvi_files <- list.files(path = "neon-nit/DP3.30018.001/2017/FullSite/D06/2017_KONZ_2/L3/Spectrometer/VegIndices/", full.names = T, recursive = T, pattern = "NDVI.tif") ndvi_rasters <- lapply(X = ndvi_files, FUN = raster) ndvi_rasters$fun <- max ndvi_rasters$na.rm <- TRUE ndvi <- do.call(what = mosaic, args = ndvi_rasters) ndvi <- mask(x = ndvi, mask = sitebnds_KONZ) ndvi <- crop(x = ndvi, y = sitebnds_KONZ) ndvi.5 <- ndvi ndvi.5[ndvi.5<=.5] <- NA writeRaster(x = ndvi.5, filename = "KONZ_mosaics/KONZ_ndvi0.5.tif") ###### full spec (Konza_flyover_ground__truth_data06082017_revis) https://urldefense.com/v3/__https://www.neonscience.org/hsi-hdf5-r__;!!PhOWcWs!kCIISsJohuXS_AzvmX9EZylKdU-H4ndnceaaHYQvLsWwOOt6Z8HZBd2FTn6n$ dir.create(path = "KONZ_mosaics/spc", showWarnings = FALSE) #cores <- 3 #c1 <- makeCluster(cores) #registerDoSNOW(c1) kzmask <- reclassify(x = chm, rcl = matrix(data = c(-Inf,Inf,1), ncol = 3)) h5file <- list.files(path = "neon-spc/DP3.30006.001/2017/FullSite/D06/2017_KONZ_2/L3/Spectrometer/Reflectance/", pattern = ".h5", full.names = TRUE, recursive = TRUE) #dir.create("./temp") bands <- 1:426 #foreach (b = bands, .packages = c('raster', 'plyr', 'rhdf5')) %dopar% { for(b in bands){ print(paste("Processing band: ", b,sep = "")) options(rasterTmpDir=paste('./temp/', b, sep="")) f.list <- rep(list(NULL),length(h5file)) for(f in h5file){ print(paste("Processing swath: ", basename(f),sep = "")) wl <- h5read(f, "/KONZ/Reflectance/Metadata/Spectral_Data")$Wavelength refl <- t(h5read(f, "/KONZ/Reflectance/Reflectance_Data", index = list(b,NULL,NULL))[1,,]) proj.neon <- sub(pattern = "UTM", replacement = "utm", x = h5read(f, "/KONZ/Reflectance/Metadata/Coordinate_System")$Proj4) ext <- h5read(f, "/KONZ/Reflectance/Metadata/Coordinate_System")$Map_Info # The 4th and 5th columns of map info signify the coordinates of the map origin, which refers to the upper-left corner of the image (xMin, yMax). ext <- unlist(strsplit(ext, split = ",\\s*")) # Extract the upper left-hand corner coordinates from mapInfo xmin <- as.numeric(ext[4]) ymax <- as.numeric(ext[5]) # Calculate the xMax and yMin values from the dimensions # xMax = left corner + (# of columns * resolution) xmax = xmin + (dim(refl)[2] * as.numeric(ext[6])) ymin = ymax - (dim(refl)[1] * as.numeric(ext[7])) ext <- extent(c(xmin, xmax, ymin, ymax)) refl <- raster(x = refl, crs = CRS(projargs = proj.neon)) extent(refl) <- ext refl[refl == -9999] <- NA f.list[[which(h5file == f)]] <- refl # H5close(f) } f.list$fun <- max f.list$na.rm <- TRUE b.mosaic <- do.call(mosaic, f.list) #clip to KZA b.mosaic <- raster::mask(crop(b.mosaic, kzmask), kzmask) fname <- paste("KONZ_mosaics/spc/", "band_", sprintf("%03d", b), ".tif", sep = "") writeRaster(x = b.mosaic, file = fname, overwrite=TRUE, datatype = "INT2S") file.remove(list.files(path = paste('./temp/', b, sep=""),full.names = T)) } #stopCluster(c1)