library(magrittr)
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:magrittr':
##
## extract
library(rgdal)
## rgdal: version: 1.2-8, (SVN revision 663)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.1, released 2017/06/23
## Path to GDAL shared files: /usr/share/gdal/2.2
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-5
library(gdalUtils)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
setwd("/home/michael/Dropbox/BGU/Meitar_Sorek_Hamer/p_09_maiac_hdf_files_to_point_layer_Alaa_Mhawish/")
# Read data from MAIAC static grid files (for analyzed tile)
grid = get_subdatasets("MAIACLatlon.h00v02.hdf")
lon = grid[grepl("latlon:lon", grid)] %>% readGDAL %>% raster
## HDF4_EOS:EOS_GRID:MAIACLatlon.h00v02.hdf:latlon:lon has GDAL driver HDF4Image
## and has 1200 rows and 1200 columns
lat = grid[grepl("latlon:lat", grid)] %>% readGDAL %>% raster
## HDF4_EOS:EOS_GRID:MAIACLatlon.h00v02.hdf:latlon:lat has GDAL driver HDF4Image
## and has 1200 rows and 1200 columns
# Add 'row' and 'col' layers
row = lon
row[] = rowFromCell(lon, 1:ncell(lon))
col = lon
col[] = colFromCell(lon, 1:ncell(lon))
# Combine to multi-band raster with 'row', 'col', 'lon' and 'lat' values
grid = stack(row, col, lon, lat)
names(grid) = c("row", "col", "lon", "lat")
plot(grid)
# Convert raster to data.frame
grid = as.data.frame(grid)
head(grid)
## row col lon lat
## 1 1 1 64.23634 32.03839
## 2 1 2 64.24648 32.04133
## 3 1 3 64.25661 32.04427
## 4 1 4 64.26675 32.04720
## 5 1 5 64.27689 32.05014
## 6 1 6 64.28703 32.05308
# Read data from MAIAC file
sds = get_subdatasets("MAIACAAOT.h00v02.20060010850.hdf")
Optical_Depth_047 = sds[grepl("grid1km:Optical_Depth_047", sds)] %>% readGDAL %>% raster
## HDF4_EOS:EOS_GRID:MAIACAAOT.h00v02.20060010850.hdf:grid1km:Optical_Depth_047 has GDAL driver HDF4Image
## and has 1200 rows and 1200 columns
Optical_Depth_055 = sds[grepl("grid1km:Optical_Depth_055", sds)] %>% readGDAL %>% raster
## HDF4_EOS:EOS_GRID:MAIACAAOT.h00v02.20060010850.hdf:grid1km:Optical_Depth_055 has GDAL driver HDF4Image
## and has 1200 rows and 1200 columns
# Note that some variables such as 'RelAZ' are on a different resolution and must be disaggregated to match -
RelAZ = sds[grepl("grid5km:RelAZ", sds)] %>% readGDAL %>% raster %>% disaggregate(5)
## HDF4_EOS:EOS_GRID:MAIACAAOT.h00v02.20060010850.hdf:grid5km:RelAZ has GDAL driver HDF4Image
## and has 240 rows and 240 columns
# More variables...
AOT_Uncertainty = sds[grepl("grid1km:AOT_Uncertainty", sds)] %>% readGDAL %>% raster
## HDF4_EOS:EOS_GRID:MAIACAAOT.h00v02.20060010850.hdf:grid1km:AOT_Uncertainty has GDAL driver HDF4Image
## and has 1200 rows and 1200 columns
AOT_QA = sds[grepl("grid1km:AOT_QA", sds)] %>% readGDAL %>% raster
## HDF4_EOS:EOS_GRID:MAIACAAOT.h00v02.20060010850.hdf:grid1km:AOT_QA has GDAL driver HDF4Image
## and has 1200 rows and 1200 columns
# Create additional 'row' and 'col' layers to join with reference grid
row = Optical_Depth_047
row[] = rowFromCell(Optical_Depth_047, 1:ncell(Optical_Depth_047))
col = Optical_Depth_047
col[] = colFromCell(Optical_Depth_047, 1:ncell(Optical_Depth_047))
# Combine to multi-band raster with 'row', 'col' and AOD values
r = stack(
row,
col,
Optical_Depth_047,
Optical_Depth_055,
RelAZ,
AOT_Uncertainty,
AOT_QA
)
names(r) = c(
"row",
"col",
"Optical_Depth_047",
"Optical_Depth_055",
"RelAZ",
"AOT_Uncertainty",
"AOT_QA"
)
plot(r)
# Convert to data.frame
r = as.data.frame(r)
head(r)
## row col Optical_Depth_047 Optical_Depth_055 RelAZ AOT_Uncertainty AOT_QA
## 1 1 1 0.246 0.182 63.67 0.0909 1
## 2 1 2 0.289 0.214 63.67 0.0974 33
## 3 1 3 0.314 0.233 63.67 0.1277 33
## 4 1 4 NA NA 63.67 0.1130 3
## 5 1 5 NA NA 63.67 0.1028 3
## 6 1 6 NA NA 65.20 0.0952 3
# Join AOD data with reference grid, to determine 'lon' and 'lat' for each AOD estimate
r = left_join(r, grid, c("row", "col"))
r = r[!is.na(r$lon) & !is.na(r$lat), ]
# Final result table
head(r)
## row col Optical_Depth_047 Optical_Depth_055 RelAZ AOT_Uncertainty AOT_QA
## 1 1 1 0.246 0.182 63.67 0.0909 1
## 2 1 2 0.289 0.214 63.67 0.0974 33
## 3 1 3 0.314 0.233 63.67 0.1277 33
## 4 1 4 NA NA 63.67 0.1130 3
## 5 1 5 NA NA 63.67 0.1028 3
## 6 1 6 NA NA 65.20 0.0952 3
## lon lat
## 1 64.23634 32.03839
## 2 64.24648 32.04133
## 3 64.25661 32.04427
## 4 64.26675 32.04720
## 5 64.27689 32.05014
## 6 64.28703 32.05308
# To point layer
coordinates(r) = ~ lon + lat
proj4string(r) = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
# Plot example
plot(r[sample(1:nrow(r), 500), ])
