[R-sig-Geo] Preprocess daily Black Marble (NVP46A2) nighttime light product

[Nikolaos Tziokas]

I downloaded NASA's Black Marble daily product (VNP46A2) which is in h5
format (the data can be dowloaded from their website
<https://ladsweb.modaps.eosdis.nasa.gov/search/order/1/VNP46A2--5000> (it's
free you just need to create and account) or from my GDrive
<https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7>).
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.

I can understand how to syntax the ifel for the Mandatory_Quality_Flag and
Snow_Flag but I can't understand how to syntax the bits for the
QF_Cloud_Mask.

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
-- 
Tziokas Nikolaos
Cartographer

Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>

	[[alternative HTML version deleted]]