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]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
