I downloaded NASA's Black Marble monthly nighttime light NTL data, VNP46A3. Ι am trying to set the projection to EPSG:4326 but the resulting product (All_Angles_Snow_Free) seems to go go way off. For example, I wanted to download NTL data for the city of Mumbai, India. After reprojecting the NTL (product 5 (All_Angles_Snow_Free) from the .h5) the result is shown in the attachement.
At the bottom if the image is a shp of Mumbai and the red circle in the top indicates where Mumbai is in the NTL image. Clearly something's not right.
I downloaded the image from here (LAADS-DAAC, Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center). The code I used to extract the NTL radiance image is:
library(terra)
wd <- "path/"
r <- rast(paste0(wd, "VNP46A3.A2018091.h25v07.001.2021125122857.h5"))
crs(r) <- "epsg:4326"
2400*(15/(60*60))
h = 25
v = 7
ext(r) = c(-180+h*10,-180+(h+1)*10, (v-2)*10,(v-1)*10)
ntl <- r[[5]]
writeRaster(ntl, paste0(wd, "ntl.tif"), overwrite = TRUE)
What's wrong with the code? From here (https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7?usp=drive_link) you can download the .h5 image if you don't want to use NASA's website. I am using R 4.3.1 and RStudio 2023.06.2+561.
NASA's Black Marble monthly data: Reprojection isn't accurate
-
- Posts: 1
- Joined: Sat Sep 09, 2023 10:58 am America/New_York
NASA's Black Marble monthly data: Reprojection isn't accurate
- Attachments
-
- difference between the reprojected image and a shp.
- PYMXr.png (5.29 KiB) Not viewed yet
Tags:
-
- User Services
- Posts: 266
- Joined: Mon Sep 30, 2019 8:33 am America/New_York
- Has thanked: 1 time
Re: NASA's Black Marble monthly data: Reprojection isn't accurate
Hi Nikolaos,
Yes, it appears that the issue lies in how you're calculating the spatial extent. This is how you calculate the bounding box using the horizontal (h) and vertical(v) time number –
WestBoundCoord = (10*HorizontalTileNumber) - 180
NorthBoundCoord = 90-(10*VerticalTileNumber)
EastBoundCoord = WestBoundCoord + 10
SouthBoundCoord = NorthBoundCoord - 10
I've noticed that you're using the R program to extract the VNP46A3 layer. We actually have a Python script that accomplishes the same task. Here is a script for you reference –
https://blackmarble.gsfc.nasa.gov/tools/OpenHDF5.py
Yes, it appears that the issue lies in how you're calculating the spatial extent. This is how you calculate the bounding box using the horizontal (h) and vertical(v) time number –
WestBoundCoord = (10*HorizontalTileNumber) - 180
NorthBoundCoord = 90-(10*VerticalTileNumber)
EastBoundCoord = WestBoundCoord + 10
SouthBoundCoord = NorthBoundCoord - 10
I've noticed that you're using the R program to extract the VNP46A3 layer. We actually have a Python script that accomplishes the same task. Here is a script for you reference –
https://blackmarble.gsfc.nasa.gov/tools/OpenHDF5.py
Regards,
LAADS User Services
To receive news from LAADS DAAC direct to your inbox, email laadsdaac-join@lists.nasa.gov with “subscribe” in the subject line.
LAADS User Services
To receive news from LAADS DAAC direct to your inbox, email laadsdaac-join@lists.nasa.gov with “subscribe” in the subject line.