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
library(maps)

setwd("/home/michael/Dropbox/BGU/Itai_Kloog/p_52_maiac_projection_data_2/")

#####################################################################################
# Method 1 - Using reference grid

# Read 'lon' and 'lat' rasters from reference grid
grid = get_subdatasets("MAIACLatlon.h04v04.hdf")
lon = grid[grepl("latlon:lon", grid)] %>% readGDAL %>% raster 
## HDF4_EOS:EOS_GRID:MAIACLatlon.h04v04.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.h04v04.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.h04v04.20160061845.hdf")
Optical_Depth = sds[grepl("grid1km:Optical_Depth_047", sds)] %>% readGDAL %>% raster 
## HDF4_EOS:EOS_GRID:MAIACAAOT.h04v04.20160061845.hdf:grid1km:Optical_Depth_047 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)

# USA layer
usa = st_as_sf(map('usa', plot = FALSE, fill = TRUE))
laea = st_crs("+proj=laea +lat_0=30 +lon_0=-95") # Lambert equal area
usa = st_transform(usa, laea)
pnt = st_transform(pnt, laea)

# Map 1
g = st_graticule(pnt)
plot(st_geometry(usa), graticule = TRUE, key.pos = NULL, axes = TRUE)
plot(st_geometry(pnt[sample(1:nrow(pnt), 1000), ]), col = "red", add = TRUE)

#####################################################################################
# Method 2 - Reading CRS info directly with GDAL

Optical_Depth = sds[grepl("grid1km:Optical_Depth_047", sds)] %>% readGDAL %>% raster 
## HDF4_EOS:EOS_GRID:MAIACAAOT.h04v04.20160061845.hdf:grid1km:Optical_Depth_047 has GDAL driver HDF4Image 
## and has 1200 rows and 1200 columns
pnt2 = rasterToPoints(Optical_Depth, spatial = TRUE)
pnt2 = st_as_sf(pnt2)
pnt2 = st_transform(pnt2, laea)

# Map 2
g = st_graticule(pnt2)
plot(st_geometry(pnt2[sample(1:nrow(pnt2), 1000), ]), graticule = TRUE, key.pos = NULL, axes = TRUE, col = "red")