Resample_example2_v02.R

library(stars)
library(starsExtra)

###############################################
# Layer A

# Read
f = "/media/qnap/Projects/P064.IL.Israel.FUSION.PM/raw/AOD_Fusion/MAIAC/MCD19A2/MCD19A2.A2011002.h20v05.006.2018048162814.hdf"
a1 = read_stars(gdal_subdatasets(f)[[1]])
f = "/media/qnap/Projects/P064.IL.Israel.FUSION.PM/raw/AOD_Fusion/MAIAC/MCD19A2/MCD19A2.A2011002.h21v05.006.2018048233245.hdf"
a2 = read_stars(gdal_subdatasets(f)[[1]])

# Average per pixel
a1 = st_apply(a1, 1:2, mean, na.rm = TRUE)
a2 = st_apply(a2, 1:2, mean, na.rm = TRUE)

# Mosaic
a = st_mosaic(a1, a2)

# Plot - layer A
names(a) = "layer A"
plot(a, col = hcl.colors(11, "Spectral"), reset = FALSE)

###############################################
# Layer B

# Read
b = read_stars("/media/qnap/Projects/P064.IL.Israel.FUSION.PM/raw/AOD_Fusion/ECMWF_AOD/macc_2011_01.ncdf", quiet = TRUE)

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

## Warning: ignoring unrecognized unit: ~

# Subset variable
b = b["aod469"]

# Subset date
times = st_get_dimension_values(b, "time")
dates = as.Date(times)
b = b[,,,which(dates == as.Date("2011-01-02"))]

# Average per pixel
b = st_apply(b, 1:2, mean)

# Set CRS
st_crs(b) = 4326

# Plot - layer B
names(b) = "layer B"
plot(b, col = hcl.colors(11, "Spectral"), reset = FALSE)
plot(b %>% st_as_sf %>% st_geometry %>% st_transform(st_crs(a)), add = TRUE)

###############################################
# Crop A according to B

# AOI
aoi = st_bbox(b)
aoi = st_as_sfc(aoi)
aoi = st_transform(aoi, st_crs(a))

# Crop
a = st_crop(a, aoi)

# Trim
a = trim(a)

# Plot - layer A cropped
names(a) = "layer A cropped"
plot(a, col = hcl.colors(11, "Spectral"), reset = FALSE)

###############################################
# Resample

# Resample B according to A
b_resampled = st_warp(b, a, method = "bilinear", use_gdal = TRUE)

# Plot - resampled layer B
names(b_resampled) = "layer B resampled"
plot(b_resampled, reset = FALSE)
plot(a, col = hcl.colors(11, "Spectral"), add = TRUE)

# Fill missing in A according to resampled B
a[[1]][is.na(a[[1]])] = b_resampled[[1]][is.na(a[[1]])]

# Plot - final result
names(a) = "final result"
plot(a, col = hcl.colors(11, "Spectral"))