code_01.R

# Tile 
# /media/qnap/Data-NoBkp/move/MAIAC_EUROPE_V12015/dataportal.nccs.nasa.gov/DataRelease/Europe/h01v01/2014/MAIACAAOT.h01v01.20140011115.hdf

# Ref
# /media/qnap/Data-NoBkp/move/MAIAC_EUROPE_V12015/MAIACLatlon.h01v01.hdf

library(magrittr)
library(rgdal)

## Loading required package: sp

## rgdal: version: 1.3-2, (SVN revision 755)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  GDAL binary built with GEOS: TRUE 
##  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.3-1

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:magrittr':
## 
##     extract

library(gdalUtils)
library(sf)

## Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2

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/Itai_Kloog/p_51_maiac_projection_data/")

# GDAL info
system("gdalinfo MAIACAAOT.h01v01.20140011115.hdf")

# Read 'lon' and 'lat' rasters from reference grid
grid = get_subdatasets("MAIACLatlon.h01v01.hdf")
lon = grid[grepl("latlon:lon", grid)] %>% readGDAL %>% raster 

## HDF4_EOS:EOS_GRID:MAIACLatlon.h01v01.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.h01v01.hdf:latlon:lat has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns

# Creare 'row' and 'col' rasters
row = lon
row[] = rowFromCell(lon, 1:ncell(lon))
col = lon
col[] = colFromCell(lon, 1:ncell(lon))

# Combine to multi-band raster
grid = stack(row, col, lon, lat)
names(grid) = c("row", "col", "lon", "lat")

# Convert to data.frame
grid = as.data.frame(grid)

# Read 'measurements' data
sds = get_subdatasets("MAIACAAOT.h01v01.20140011115.hdf")
Optical_Depth = sds[grepl("grid1km:Optical_Depth_Land", sds)] %>% readGDAL %>% raster 

## HDF4_EOS:EOS_GRID:MAIACAAOT.h01v01.20140011115.hdf:grid1km:Optical_Depth_Land has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns

row = Optical_Depth
row[] = rowFromCell(Optical_Depth, 1:ncell(Optical_Depth))
col = Optical_Depth
col[] = colFromCell(Optical_Depth, 1:ncell(Optical_Depth))
r = stack(
  row, 
  col, 
  Optical_Depth
)
names(r) = c(
  "row", 
  "col", 
  "Optical_Depth"
)
r = as.data.frame(r)

# Join measurements with reference grid
r = left_join(r, grid, c("row", "col"))

# To layer
pnt = st_as_sf(r, coords = c("lon", "lat"), crs = 4326)

# Map
g = st_graticule(pnt)
plot(pnt[sample(1:nrow(pnt), 10000), ], graticule = TRUE, key.pos = NULL, axes = TRUE)