2

I am now working with VIIRS/NPP Active Fires by using python gdal. But I cannot read the data inside the files.
Filename = "NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.hdf"
Here is some information that I can get by using gdal.

gdal.RasterCount = 0

And gdalinfo:

gdalinfo NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.hdf
Driver: HDF4/Hierarchical Data Format Release 4
Files: NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.hdf
Size is 512, 512
Coordinate System is `'
Metadata:
  AlgorithmType=OPS
  Beginning_Time_IET=[1.7056441e+15]
  BeginningTime=060027.600000Z
  DayNightFlag=Day
  EastBoundingCoord=123.866
  Ending_Time_IET=[1.7056444e+15]
  EndingTime=060609.000000Z
  EndTime=2012-01-19 06:06:09.000
  HDFEOSVersion=HDFEOS_V2.17
InputPointer=NPP_GRCMAE_L1.A2012019.0555.P1_03110.2014057115525.hdf,NPP_GRCMAE_L1.A2012019.0600.P1_03110.2014057115623.hdf,NPP_GRCMAE_L1.A2012019.0605.P1_03110.2014057115525.hdf
  InstrumentShortname=VIIRS
  LocalGranuleID=NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.hdf
  LongName=VIIRS/NPP Active Fires 5-Min L2 Swath ARP 750m
  LPEATE_AlgorithmVersion=NPP_PRVAF 1.5.07.01
  LUTs_used=VIIRS-AF-EDR-AC-Int_v1.5.06.02_LP
  NorthBoundingCoord=27.8258
  Number_Fire_Pixels=256
  NumSCEA_RDR_TimeSegments=[18]
  NumSci_RDR_TimeSegments=[4]
  PGE_EndTime=2012-01-19 06:05:00.000
  PGE_Name=PGE330
  PGE_StartTime=2012-01-19 06:00:00.000
  PGEVersion=P2.3.0
  Platform_Short_Name=NPP
  ProcessingEnvironment=Linux minion5609 2.6.18-371.1.2.el5 #1 SMP Tue Oct 22 12:51:53 EDT 2013 x86_64 x86_64 x86_64 GNU/Linux
  ProcessVersion=P1_03110
  ProductionTime=2014-02-26 12:59:56.000
  ProxyDataType=Operational Data
  Resolution=Imagery
  SatelliteInstrument=NPP_OPS
  ShortName=NPP_AVAF_L2
  SouthBoundingCoord=3.93342
  StartTime=2012-01-19 06:00:27.600
  Unagg_DayNightFlag=TS 0: Day; TS 1: Day; TS 2: Day; TS 3: Day
  WestBoundingCoord=90.5542
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

This is what I can see from HDFViewer on same file: enter image description here

I know that the file carries information of active fire points, but I cannot read the data inside into an array.
From gdalinfo above, I can see Number_Fire_Pixels=256, but how can I get lat,lon of these points?

How can I read the data by using python?

UPDATE: Here's the link to the file that I'm working with.

  • Does gdalinfo not list any subdatasets? – mdsumner Apr 8 '15 at 4:17
  • As you can see in gdalinfo, it says Number_Fire_Pixels=256. How can I get lat, lon of these points? – Jackie Apr 8 '15 at 4:27
  • 2
    If there's no listing of subdatasets (and otherwise no primary dataset as it seems here), you can't get anything with this build of GDAL. This is essentially a "header-only" summary I'm afraid. You might have to use the HDF4 library more directly, or dump out stuff with command line utilites. Can you point to an example file? – mdsumner Apr 8 '15 at 4:34
  • GDAL cannot read these 1D arrays (or 2D with a degenerate dimension, whatever shakes your tree). You'll need the hdf library in Python (no experience here). – mdsumner Apr 8 '15 at 8:25
  • Please suggest me one library that work with python 2.7, it would be even better with instruction. Thanks in advance! – Jackie Apr 8 '15 at 13:51
3

Your data are stored as tables rather than gridded (raster) data which could be interpreted by GDAL.

It might be easier in the end to work in HDF5 rather than HDF4. Given you're on a Windows box it's easy to download and install the h4toh5 tools from the HDF group which can be used from the command line with (using your example file):

h4toh5convert NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.hdf NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.h5

Once you've got the file in HDF5 (you should not lose any data here as HDF4 attributes translate directly to HDF5) you've got a number of options in Python for working with the data. Your best bet is to install PyTables and Pandas and use the inbuilt HDFStore object to read in the data, which might look something like:

import pandas as pd

path = "NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.h5"
store = pd.io.pytables.HDFStore(path)

print store

Printing the store gives you a list of the series and frames stored in the dataset:

<class 'pandas.io.pytables.HDFStore'>
File path: /TEMP/NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.h5
/ActiveFires_ARR/Data Fields/ColIndex                  frame_table [0.0.0] (typ->generic,nrows->256,ncols->1,indexers->[index],dc->[ColIndex])      
/ActiveFires_ARR/Data Fields/Latitude                  frame_table [0.0.0] (typ->generic,nrows->256,ncols->1,indexers->[index],dc->[Latitude])      
/ActiveFires_ARR/Data Fields/Longitude                 frame_table [0.0.0] (typ->generic,nrows->256,ncols->1,indexers->[index],dc->[Longitude])     
/ActiveFires_ARR/Data Fields/QF1_VIIRSAFARP            frame_table [0.0.0] (typ->generic,nrows->256,ncols->1,indexers->[index],dc->[QF1_VIIRSAFARP])
/ActiveFires_ARR/Data Fields/QF2_VIIRSAFARP            frame_table [0.0.0] (typ->generic,nrows->256,ncols->1,indexers->[index],dc->[QF2_VIIRSAFARP])
/ActiveFires_ARR/Data Fields/QF3_VIIRSAFARP            frame_table [0.0.0] (typ->generic,nrows->256,ncols->1,indexers->[index],dc->[QF3_VIIRSAFARP])
/ActiveFires_ARR/Data Fields/QF4_VIIRSAFARP            frame_table [0.0.0] (typ->generic,nrows->256,ncols->1,indexers->[index],dc->[QF4_VIIRSAFARP])
/ActiveFires_ARR/Data Fields/RowIndex                  frame_table [0.0.0] (typ->generic,nrows->256,ncols->1,indexers->[index],dc->[RowIndex])  

Which you can then select out the data using the name as a key:

latitude = store["/ActiveFires_ARR/Data Fields/Latitude"]
longitude = store["/ActiveFires_ARR/Data Fields/Longitude"]
qf1 = store["/ActiveFires_ARR/Data Fields/QF1_VIIRSAFARP"]
qf2 = store["/ActiveFires_ARR/Data Fields/QF2_VIIRSAFARP"]
qf3 = store["/ActiveFires_ARR/Data Fields/QF3_VIIRSAFARP"]
qf4 = store["/ActiveFires_ARR/Data Fields/QF4_VIIRSAFARP"]

And then join into a single pandas dataframe object:

df = latitude.join([longitude, qf1, qf2, qf3, qf4])

And from there you've got a complete dataframe object where you can do what you like, eg:

import seaborn
df.plot(kind="scatter",
    x="Longitude", y="Latitude",
    s=df["QF4_VIIRSAFARP"] / df["QF4_VIIRSAFARP"].max() * 200,
    alpha=0.5
)

Plot of QF4 data by latitude and longitude

And one last note, the store object can be used as a context manager so it's closed automatically:

with HDFStore(path) as store:
    # do something...
  • I updated the question with link to my file. Thanks! – Jackie Apr 9 '15 at 2:58
  • @Jackie Thanks for that! I've updated the answer with more information based on your sample – om_henners Apr 9 '15 at 23:22
  • @om_henners Thanks for your solution. One more thing, how can I pick value at particular index?. For example: print longitude[index] – Jackie Apr 10 '15 at 8:01
  • To be clear, I want to convert longitude into an array, where I can work with the data. By the way, .iloc[index] can get the selected index. – Jackie Apr 10 '15 at 8:17
  • 1
    @Jackie Because pandas is build on top of numpy arrays if you select a column (or columns) from the dataframe you can cast it to a numpy array. For example numpy.array(df["Longitude"]) or numpy.array(df[["Latitude", "Longitude"]]) – om_henners Apr 10 '15 at 11:54
6

Think hdf file as a folder. You want to open the file INSIDE the folder.

import gdal
hdf_file = gdal.Open("3B43.20140501.7.HDF") # 3b43 rainfall dataset

subDatasets = hdf_file.GetSubDatasets()

subDatasets 
>>> [('HDF4_SDS:UNKNOWN:"3B43.20140501.7.HDF":0', '[1440x400] precipitation (32-bit floating-point)'), ('HDF4_SDS:UNKNOWN:"3B43.20140501.7.HDF":1', '[1440x400] relativeError (32-bit floating-point)'), ('HDF4_SDS:UNKNOWN:"3B43.20140501.7.HDF":2', '[1440x400] gaugeRelativeWeighting (8-bit integer)')]

# Open precipitation
# prcp = gdal.Open('HDF4_SDS:UNKNOWN:"3B43.20140501.7.HDF":0')
# or the following shortcut:
prcp = gdal.Open(subDatasets[0][0])
prcp.ReadAsArray()

array([[ 0.055     ,  0.07040323,  0.04701613, ...,  0.06721774,
         0.07008065,  0.07181452],
       [ 0.06096774,  0.07983872,  0.09064516, ...,  0.07157258,
         0.07733872,  0.07399193],
       [ 0.0703629 ,  0.08100805,  0.09028225, ...,  0.07931452,
         0.08270162,  0.08221775],
       ..., 
       [ 0.04266129,  0.02157258,  0.03274193, ...,  0.08129031,
         0.07431452,  0.07338709],
       [ 0.0278629 ,  0.02370968,  0.04048387, ...,  0.07133064,
         0.07189515,  0.07112902],
       [ 0.03225806,  0.03040322,  0.03907258, ...,  0.0716129 ,
         0.07233871,  0.07262097]], dtype=float32)

For a more in-depth turorial you can read http://jgomezdans.github.io/gdal_notes/ipython.html

  • gdalinfo would have listed the subdatasets if they are there, this is not that kind of file afaics - still worth seeing what python returns from GetSubDatasets() – mdsumner Apr 8 '15 at 6:52
  • It returns nothing when I try this print hdf.GetSubDatasets() – Jackie Apr 8 '15 at 6:57
  • 2
    a sample would help to clarify if you provide one – nickves Apr 8 '15 at 7:52
  • @nickves I added my file. Please have a look. Thank you! – Jackie Apr 9 '15 at 2:58
2

My script uses a NDVI (no corrected by scale factor) sub dataset of modis product, for getting the coordinates (sinusoidal projection) for a value of 256 (equivalent to Number_Fire_Pixels = 256):

from osgeo import gdal
import struct

nameraster = "MOD13Q1.A2005193.h10v08.005.2008215173619.hdf"

hdf_file = gdal.Open(nameraster)

subDatasets = hdf_file.GetSubDatasets()

dataset = gdal.Open(subDatasets[0][0])
geotransform = dataset.GetGeoTransform()
band = dataset.GetRasterBand(1)

fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 
            'Int32':'i', 'Float32':'f', 'Float64':'d'}

print "rows = %d columns = %d" % (band.YSize, band.XSize)

BandType = gdal.GetDataTypeName(band.DataType)

print "Data type = ", BandType

print "Executing with %s" % nameraster

print "test_value = 256"

X = geotransform[0] #x coordinate
Y = geotransform[3] #y coordinate

for y in range(band.YSize):

    scanline = band.ReadRaster(0, 
                               y, 
                               band.XSize, 
                               1, 
                               band.XSize, 
                               1, 
                               band.DataType)

    values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)

    for value in values:

        if(value == 256):       
            print "%.4f %.4f %.2f" % (X, Y, value)
        X += geotransform[1] #x pixel size
    X = geotransform[0]
    Y += geotransform[5] #y pixel size

dataset = None

When I run the script at the bash shell I get:

rows = 4800 columns = 4800
Data type =  Int16
Executing with MOD13Q1.A2005193.h10v08.005.2008215173619.hdf
test_value = 256
-8188820.6083 1048476.6775 256.00
-8188820.6083 1048245.0211 256.00
-7870988.0847 1031797.4197 256.00
-8567115.4413 1028554.2307 256.00
-8526343.9223 981296.3336 256.00
-8739236.1155 975968.2374 256.00
-8180249.3230 972493.3920 256.00
-8690819.9366 957899.0414 256.00
-8771668.0057 954192.5397 256.00
-8703329.3800 915274.2715 256.00
-8585416.2936 915042.6151 256.00
-7988669.5147 722999.4941 256.00
-7986584.6075 717208.0852 256.00
-7994924.2364 699833.8583 256.00
-7992839.3292 697748.9511 256.00
-7992607.6728 697748.9511 256.00
-7993534.2983 696822.3257 256.00
-8177932.7594 553658.6963 256.00
-8177701.1031 553658.6963 256.00
-7904114.9440 532114.6549 256.00
-8148975.7146 505474.1737 256.00
-8148744.0583 505242.5174 256.00
-8595145.8607 460532.8402 256.00
-8595145.8607 460301.1839 256.00
-8410979.0558 408873.4723 256.00
-8410747.3995 408873.4723 256.00
-8410979.0558 408641.8160 256.00
-8410747.3995 408641.8160 256.00
-8408199.1795 407715.1905 256.00
-8429743.2209 372040.1114 256.00
-8441094.3824 325477.1834 256.00
-8442484.3206 325013.8706 256.00
-8658156.3901 293045.2932 256.00
-8764486.6586 153819.8219 256.00
-8764255.0022 153819.8219 256.00
-8764486.6586 153588.1655 256.00
-8764255.0022 153588.1655 256.00
-8317621.5435 146638.4748 256.00
-8377620.5403 85481.1962 256.00

Editing Note:

I have in my GNU/Linux Debian a hdf visor (HDFView). I opened the NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.hdf file and I accessed the data with a simple copy/paste in a spreadsheet. Afterward, I saved them as a *.cvs file and loaded the file at QGIS. This was the result:

enter image description here

This is the image of NPP_AVAF_L2.A2012019.0600.P1_03110.2014057125956.hdf file opened in HDFView. Icons of each dataset seems to point out that they are databases.

enter image description here

Below image is of the modis product used with my script. In this case, icons of each dataset point out that they are grids.

enter image description here

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.