## Install packages
## install.packages("gdalUtils")
## install.packages("raster")
## install.packages("rgdal")
## install.packages("rgeos")
## install.packages("plyr")
## install.packages("magrittr")
## install.packages("sp")
## install.packages("ncdf4")
## install.packages("reshape2")
## install.packages("foreach")
## install.packages("parallel")
## install.packages("doParallel")
## install.packages("mapview")
## install.packages("velox")
## install.packages("sf")
###########################################################
## Part a.1- working with hdf files on Windows ############
###########################################################
suppressMessages(library(gdalUtils))
suppressMessages(library(raster))
suppressMessages(library(rgdal))
suppressMessages(library(rgeos))
## Warning: package 'rgeos' was built under R version 3.4.4
suppressMessages(library(plyr))
suppressMessages(library(magrittr))
suppressMessages(library(sp))
suppressMessages(library(mapview))
# If you do not have it installed, you will probably get an error message and won't be able to read hdf files
# gdal_setInstallation()
# x <- getOption("gdalUtils_gdalPath")
# x[[1]]$path ## Look where GDAL is installed
# df <- x[[1]]$drivers # Table of all installed drivers
# Set your working directory
setwd("D:/MAIAC_MCD19")
###########################################################################
# For more information about the MAIAC AOD product and the different layers of each subdatasets, go to:
# https://lpdaac.usgs.gov/node/1265
# https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/maiac/MCD19A2/
################# Test on one HDF file ####################################
# Get a list of sds names
# sds <- get_subdatasets("MCD19A2.A2018183.h13v11.006.2018186225614.hdf") # Returns subdataset names
#
# # R wrapper for gdal_translate - converts from hdf to tiff file
# gdal_translate(sds[1], "test.tif")
#
# # Load the multi-layer raster into R
# r = brick("test.tif")
# plot(r)
# mapview(r)
###########################################################################
# This function is useful if you are running the code on Windows
read_hdf = function(file, n) {
sds = get_subdatasets(file)
f = tempfile(fileext = ".tif") # the tiff file is saved as a temporary file
gdal_translate(sds[n], f)
brick(f) # the ouptut is a multi-layer raster file
}
# Define input Directories
aod_dir = "D:/MAIAC_MCD19"
#####################
for(year in 2018) {
# Read HDF files list from AOD directory
setwd(file.path(aod_dir, year))
files = list.files(pattern = "MCD19A2.*\\.hdf$", recursive = TRUE) # note that MAIACTAOT is for TERRA data and MAIACAAOT is for AQUA data
result = list()
for(f in files) {
# Read data
sds = get_subdatasets(f)
# Choose which subdatasets you want to retrieve from the hdf file
Optical_Depth_047 = read_hdf(f, grep("grid1km:Optical_Depth_047", sds))
Optical_Depth_055 = read_hdf(f, grep("grid1km:Optical_Depth_055", sds))
AOT_Uncertainty = read_hdf(f, grep("grid1km:AOD_Uncertainty", sds))
AOT_QA = read_hdf(f, grep("grid1km:AOD_QA", sds))
RelAZ=read_hdf(f, grep("grid5km:RelAZ", sds))
RelAZ=disaggregate(RelAZ, fact = 5)
# Create a different name for each layer
names(Optical_Depth_047) = paste0("Optical_Depth_047.", letters[1:nlayers(Optical_Depth_047)])
names(Optical_Depth_055) = paste0("Optical_Depth_055.", letters[1:nlayers(Optical_Depth_055)])
names(AOT_Uncertainty) = paste0("AOT_Uncertainty.", letters[1:nlayers(AOT_Uncertainty)])
names(AOT_QA) = paste0("AOT_QA.", letters[1:nlayers(AOT_QA)])
names(RelAZ) = paste0("RelAZ.", letters[1:nlayers(RelAZ)])
# Stack all the raster together
r = stack(Optical_Depth_047, Optical_Depth_055, AOT_Uncertainty, AOT_QA, RelAZ)
r = as.data.frame(r, xy=TRUE)
# Add filename
r$date =
f %>%
strsplit("\\.") %>%
sapply("[", 2) %>%
substr(2, 8) %>%
as.Date(format = "%Y%j")
# Combine results
result[[f]] = r
}
result = do.call(plyr::rbind.fill, result)
setwd(file.path(aod_dir))
saveRDS(result, sprintf("MAIACAOD_Sao_Paolo_%s.rds", year)) # Save each year separatly
}
head(result)
## x y Optical_Depth_047.a Optical_Depth_047.b
## 1 -5559289 -2224364 0.095 NA
## 2 -5558363 -2224364 0.107 NA
## 3 -5557436 -2224364 0.087 NA
## 4 -5556509 -2224364 0.092 NA
## 5 -5555583 -2224364 0.071 NA
## 6 -5554656 -2224364 0.076 NA
## Optical_Depth_047.c Optical_Depth_055.a Optical_Depth_055.b
## 1 NA 0.065 NA
## 2 NA 0.073 NA
## 3 NA 0.060 NA
## 4 NA 0.063 NA
## 5 NA 0.048 NA
## 6 NA 0.052 NA
## Optical_Depth_055.c AOT_Uncertainty.a AOT_Uncertainty.b
## 1 NA 0.0195 NA
## 2 NA 0.0199 NA
## 3 NA 0.0193 NA
## 4 NA 0.0200 NA
## 5 NA 0.0202 NA
## 6 NA 0.0192 NA
## AOT_Uncertainty.c AOT_QA.a AOT_QA.b AOT_QA.c RelAZ.a RelAZ.b RelAZ.c
## 1 NA 1 1280 1280 116.88 NA NA
## 2 NA 1 1280 1280 116.88 NA NA
## 3 NA 1 1280 1280 116.88 NA NA
## 4 NA 1 1280 1280 116.88 NA NA
## 5 NA 1 1280 1280 116.88 NA NA
## 6 NA 1 1280 1280 111.93 NA NA
## date
## 1 2018-07-01
## 2 2018-07-01
## 3 2018-07-01
## 4 2018-07-01
## 5 2018-07-01
## 6 2018-07-01
tail(result)
## x y Optical_Depth_047.a Optical_Depth_047.b
## 2879995 -4452899 -3335388 NA NA
## 2879996 -4451972 -3335388 NA NA
## 2879997 -4451045 -3335388 NA NA
## 2879998 -4450119 -3335388 NA NA
## 2879999 -4449192 -3335388 NA NA
## 2880000 -4448265 -3335388 NA NA
## Optical_Depth_047.c Optical_Depth_055.a Optical_Depth_055.b
## 2879995 NA NA NA
## 2879996 NA NA NA
## 2879997 NA NA NA
## 2879998 NA NA NA
## 2879999 NA NA NA
## 2880000 NA NA NA
## Optical_Depth_055.c AOT_Uncertainty.a AOT_Uncertainty.b
## 2879995 NA NA NA
## 2879996 NA NA NA
## 2879997 NA NA NA
## 2879998 NA NA NA
## 2879999 NA NA NA
## 2880000 NA NA NA
## AOT_Uncertainty.c AOT_QA.a AOT_QA.b AOT_QA.c RelAZ.a RelAZ.b
## 2879995 NA 11 NA 11 113.44 NA
## 2879996 NA 11 NA 11 113.28 NA
## 2879997 NA 11 NA 11 113.28 NA
## 2879998 NA 11 NA 11 113.28 NA
## 2879999 NA 11 NA 11 113.28 NA
## 2880000 NA 11 NA 11 113.28 NA
## RelAZ.c date
## 2879995 112.12 2018-07-02
## 2879996 112.36 2018-07-02
## 2879997 112.36 2018-07-02
## 2879998 112.36 2018-07-02
## 2879999 112.36 2018-07-02
## 2880000 112.36 2018-07-02