1

I'm trying to analyze nighttime lights of a specific area with VNP46A1 data. I got my data from this website, but I could find only few tutorial about this product, and one of them was VIIRS_SR_Tutorial. When I listed VNP46A1 Metadata, I found my first challenge.

from h5py import File
data_path = 'path-to-my-h5-file'
file         = File(data_path)
fileMetadata = file['HDFEOS INFORMATION/StructMetadata.0'][()].split()
fileMetadata = [m.decode('utf-8') for m in fileMetadata]
print(fileMetadata[:11])

>>>'GROUP=SwathStructure',
   'END_GROUP=SwathStructure',
   'GROUP=GridStructure',
   'GROUP=GRID_1',
   'GridName="VNP_Grid_DNB"',
   'XDim=2400',
   'YDim=2400',
   'UpperLeftPointMtrs=(20000000.000000,50000000.000000)',
   'LowerRightMtrs=(30000000.000000,40000000.000000)',
   'Projection=HE5_GCTP_GEO',
   'GridOrigin=HE5_HDFE_GD_UL'

In this page, It says "VIIRS/NPP Daily Gridded Day Night Band 500m Linear Lat Lon Grid Night" but in Production Information section it says the data has 750m spatial resolution. However, when I checked this value by dividing sides by pixel counts I found this

(30000000.000000 - 20000000.000000) / 2400
>>> 4166.666666666667

(50000000.000000 - 40000000.000000) / 2400
>>> 4166.666666666667

what is the real spatial resolution of my data ?

Also another challenge for me is was projection types. In the tutorial, the data was clearly projected Sinusoidal coordinate system,that can be seen in metadata of VNP09GA

Projection=HE5_GCTP_SNSOID

and the author defined it like this

prj = 'PROJCS["unnamed",\
GEOGCS["Unknown datum based upon the custom spheroid", \
DATUM["Not specified (based on custom spheroid)", \
SPHEROID["Custom spheroid",6371007.181,0]], \
PRIMEM["Greenwich",0],\
UNIT["degree",0.0174532925199433]],\
PROJECTION["Sinusoidal"], \
PARAMETER["longitude_of_center",0], \
PARAMETER["false_easting",0], \
PARAMETER["false_northing",0], \
UNIT["Meter",1]]'

But in my data, it says the projection of VNP46A1 is

Projection=HE5_GCTP_GEO

In this Manual, section 8.3.1 (GCTP Projection Codes), I found out what HE5_GCTP_GEO stands for, which is Geographic. But still don't know how to define it.

If the the projection, used to define my satellite data, was Geographic, are my projection coordinates and latitude-longitude values are the same? In what coordinate system my corner points are defined?

UpperLeftPointMtrs=(20000000.000000,50000000.000000)
LowerRightMtrs=(30000000.000000,40000000.000000)

I need to georeference my data but with all these unknowns I cannot do that

  • I just reached out to their support team yesterday regarding this exact issue, also referred them to this post. The metadata is insufficient, and based on what you're displaying potentially inaccurate. I'll get back to you based on what I learn. – Jake Shermeyer Aug 2 '19 at 20:54
  • This does not provide an answer to the question. Once you have sufficient reputation you will be able to comment on any post; instead, provide answers that don't require clarification from the asker. - From Review – Marcelo Villa Aug 2 '19 at 21:54
  • @JakeShermeyer it would be very helpful if you do that, thanks – gokdumano Aug 4 '19 at 21:46
  • There also the following metadata inside h5 file :NorthBoundingCoord, SouthBoundingCoord, EastBoundingCoord, WestBoundingCoord – Leonidas Sep 15 '19 at 20:16
0

I managed to figure out what was going on here after several emails with their support team. The projection is indeed a simple geographic lat/long projection with units in decimal degrees. WGS84 is the datum.

The extent labels are actually wrong. The units are not in meters they are thankfully in degrees, minutes, seconds.

'UpperLeftPointMtrs=(20000000.000000,50000000.000000)', 'LowerRightMtrs=(30000000.000000,40000000.000000)',

Conversion for this would be: 'UpperLeftPointDMS=(20 degrees 000 minutes 000.000000 seconds,50 D 000 M 000.000000 S)', 'LowerRightDMS=(30 degrees 000 minutes 000.000000 seconds,40 D 000 M 000.000000 S)'.

If you then do the math you did where you divide by 2400 and a rough conversion to meters, each cell should be close to 500m.

  • I also happen to learn that these values are not in meter but in degree. Wish there was a way to be sure about the datum tho. I'll keep looking to find it, thanks a lot. – gokdumano Aug 6 '19 at 16:56
  • @ShayzeOn1 just confirmed the datum is WGS84 – Jake Shermeyer Aug 7 '19 at 15:57
  • it is a relief to finally know that, I am sending countless e-mail but none of them sent a reply. Where did you get this info? This is my graduation project and I'm sure the professor will ask me this question. – gokdumano Aug 7 '19 at 17:46
  • @ShayzOn1, I emailed with: support-laads@earthdata.nasa.gov Also if you're using this website: ladsweb.modaps.eosdis.nasa.gov to order data, note that the grid they show is inaccurate for this product. You should pick a specific lat/long, country, or draw a box to complete an order over your AOI. Good luck with the project. The fact that this information isn't easily accessible publicly is shockingly lazy data management. – Jake Shermeyer Aug 7 '19 at 18:35
0

Follow on answer to this. Indeed the VNP46A1 suite is poorly projected and metadata is opaque and incorrect BUT NASA has provided an HDF to GeoTIFF converter tool tailor-made for rectifying this issue with Black Marble data. You can find the python code for converting HDF4 or HDF5 to GeoTIFF here: https://blackmarble.gsfc.nasa.gov/Tools.html

Just click on the "HERE" links for downloading source code for whichever converter you want (HDF 4 vs. HDF5) and follow along the steps. You may need to make slight alterations to customize to your home directory and file structure but after ~5 minutes of playing with this I was able to correctly geo-reference and project my VNP46A1 .h5 file downloaded from LAADS.

A couple things to look out for:

1) The rasterFilePre = rasterFiles[0][:-4] line has default 0 and :-4 index arguments but you will need to change the 0 to correspond to the correct file index in your folder and :-4 assumes that file extensions are three characters (plus a period). However, for my .h5 file types I had to change this to :-3

2) You will need to change the first [0] index in hdflayer.GetSubDatasets()[0][0] to correspond to the layer in the HDF file you want to extract (the 5th layer for most of the VNP46A1.h5 files corresponds to the "DNB_At_Sensor_Radiance_500m" layer FYI, so you'd reset the [0] to [4]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.