I downloaded NASA's Black Marble daily product (VNP46A2) which is in h5 format (the data can be dowloaded from their website (it's free you just need to create and account) or from here). One needs to preprocess the data using the Scientific Data Sets (SDS) included in the h5 file. Based on the User Guide, these are the following parameters I need to account for:
Table 4, page 14, Value of QF_Cloud_Mask in the VNP46A1/VJ146A1 product:
| Bit | Flag description key | Interpretation |
|---|---|---|
| 0 | Day/night | 0 = Night |
| 4-5 | Cloud Mask Quality | 11 = High |
| 6-7 | Cloud Detection Results & Confidence Indicator | 00 = Confident Clear |
| 8 | Shadow Detected | 0 = No |
| 9 | Cirrus Detection (IR) (BTM15 – BTM16) | 0 = No Cloud |
| 10 | Snow/ Ice Surface | 0 = No Snow/Ice |
Table 7, page 17, Values of the Mandatory_Quality_Flag in VNP46A2/VJ146A2 product:
| Value | Retrieval quality | Algorithm instance |
|---|---|---|
| 00 | High-quality | Main algorithm (Persistent nighttime lights) |
Table 8, page 17, Values of the Snow_Flag in VNP46A2/VJ146A2 product:
| Flag description key | Value | Interpretation |
|---|---|---|
| Snow/ Ice Surface | 00 | No Snow/Ice |
In the above tables I included the bit values I want to use to preprocess the NTL product, called Gap_Filled_DNB_BRDFCorrected_NTL.
I am using R's terra package to preprocess the product. So far what I have managed to do is:
library(terra)
wd <- "path/"
r <- rast(paste0(wd, "VNP46A2.A2018038.h28v07.001.2020333204506.h5"))
crs(r) <- "epsg:4326"
# dimensions
2400*(15/(60*60))
h = 28
v = 7
ext(r) = c(-180+h*10,-180+(h+1)*10, (8-v)*10,(8-v+1)*10) # up to this point the code works well
# the tif images inside the h5 file (for the ifel function below)
ntl <- r[[3]] # this is the Gap_Filled_DNB_BRDFCorrected_NTL
latest_high_quality_retrieval <- r[[4]]
mandatory_quality_flag <- r[[5]]
qf_cloud_mask <- r[[6]]
snow_mask <- r[[7]]
# here is the issue!!!
result <- ifel(r[[4]] > 0 & r[[5]] == 00 & r[[6]] == 1 & r[[7]] == 00, r[[3]], NA)
# scale factor based on the User Guide table 6, page 16
result1 <- result * 0.1
writeRaster(result1, paste0(wd, "ntl.tif"), overwrite = TRUE)
The writeRaster function returns an empty raster with null values.
Could you help me syntax the ifel function properly using the bits from the tables? In a very abstract sense, the ifel statement should say:
If
snow_flag is 00 AND
Mandatory_Quality_Flag is 00 AND
the bit 0 from the QF_Cloud_Mask is 0 AND
the bit 4-5 from the QF_Cloud_Mask is 11 AND
the bit 6-7 from the QF_Cloud_Mask is 0 AND
the bit 8 from the QF_Cloud_Mask is 0 AND
the bit 9 from the QF_Cloud_Mask is 0 AND
the bit 10 from the QF_Cloud_Mask is 0 THEN
keep the values of the Gap_Filled_DNB_BRDFCorrected_NTL ELSE
NA.