来自hdf文件python的lat，lon信息

Question

我有一个hdf文件，想从中提取数据。由于某些原因，我无法提取纬度和经度值：

我尝试的代码是：

from pyhdf import SD
hdf = SD.SD('MOD10C2.A2001033.006.2016092173057.hdf')
data = hdf.select('Eight_Day_CMG_Snow_Cover')
lat = (hdf.select('Latitude'))[:]

它给了我一个错误：

^{pr2}$

我试过：

lat = (hdf.select('Lat'))[:]

还是没用！在

数据可以在link中找到

非常感谢任何帮助！在

数据格式如下：

我得到的错误是：

---------------------------------------------------------------------------
HDF4Error                                 Traceback (most recent call last)
~/anaconda3/lib/python3.6/site-packages/pyhdf/SD.py in select(self, name_or_index)
   1635             try:
-> 1636                 idx = self.nametoindex(name_or_index)
   1637             except HDF4Error:

~/anaconda3/lib/python3.6/site-packages/pyhdf/SD.py in nametoindex(self, sds_name)
   1528         sds_idx = _C.SDnametoindex(self._id, sds_name)
-> 1529         _checkErr('nametoindex', sds_idx, 'non existent SDS')
   1530         return sds_idx

~/anaconda3/lib/python3.6/site-packages/pyhdf/error.py in _checkErr(procName, val, msg)
     22             err = "%s : %s" % (procName, msg)
---> 23         raise HDF4Error(err)

HDF4Error: nametoindex : non existent SDS

During handling of the above exception, another exception occurred:

HDF4Error                                 Traceback (most recent call last)
<ipython-input-11-21e6a4fdf8eb> in <module>()
----> 1 hdf.select('Lat')

~/anaconda3/lib/python3.6/site-packages/pyhdf/SD.py in select(self, name_or_index)
   1636                 idx = self.nametoindex(name_or_index)
   1637             except HDF4Error:
-> 1638                 raise HDF4Error("select: non-existent dataset")
   1639         id = _C.SDselect(self._id, idx)
   1640         _checkErr('select', id, "cannot execute")

HDF4Error: select: non-existent dataset

Answer 1

正常的经纬度信息不是hdf文件的科学模式，这是主要原因是lat = (hdf.select('Lat'))[:]不像其他变量那样工作。使用以下函数，您可以提取hdf文件中任何类型的变量存储

from pyhdf.HDF import *
from pyhdf.V   import *
from pyhdf.VS  import *
from pyhdf.SD  import *

def HDFread(filename, variable, Class=None):
    """
    Extract the data for non scientific data in V mode of hdf file
    """
    hdf = HDF(filename, HC.READ)

    # Initialize the SD, V and VS interfaces on the file.
    sd = SD(filename)
    vs = hdf.vstart()
    v  = hdf.vgstart()

    # Encontrar el puto id de las Geolocation Fields
    if Class == None:
        ref = v.findclass('SWATH Vgroup')
    else:
        ref = v.findclass(Class)

    # Open all data of the class
    vg = v.attach(ref)
    # All fields in the class
    members = vg.tagrefs()

    nrecs = []
    names = []
    for tag, ref in members:
        # Vdata tag
        vd = vs.attach(ref)
        # nrecs, intmode, fields, size, name = vd.inquire()
        nrecs.append(vd.inquire()[0]) # number of records of the Vdata
        names.append(vd.inquire()[-1])# name of the Vdata
        vd.detach()

    idx = names.index(variable)
    var = vs.attach(members[idx][1])
    V   = var.read(nrecs[idx])
    var.detach()
    # Terminate V, VS and SD interfaces.
    v.end()
    vs.end()
    sd.end()
    # Close HDF file.
    hdf.close()

    return np.array(V).ravel()

如果不知道确切的变量名，可以使用下面的程序来尝试whit pyhdf.V，该程序显示了其中包含的vgroups的内容 任何HDF文件。在

^{pr2}$

Answer 2

您使用的数据文件是MODIS 3级产品。所有3级产品都被插值到一些规则网格上。对于MOD10C2，网格就是所谓的气候模拟网格（CMG）。此网格以0.05度的间隔有规律地间隔。帕诺普知道这一点。在

CMG是柱面投影中的规则矩形网格。我们可以利用这些信息来定位数据。考虑下面的例子。在

from pyhdf import SD
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import LinearSegmentedColormap
hdf = SD.SD('MOD10C2.A2001033.006.2016092173057.hdf')
data = hdf.select('Eight_Day_CMG_Snow_Cover')
snowcover=np.array(data[:,:],np.float)
snowcover[np.where(snowcover==255)]=np.nan
m = Basemap(projection='cyl', resolution = 'l',
    llcrnrlat=-90, urcrnrlat=90,llcrnrlon=-180,urcrnrlon=180)
cdict = {'red' : [(0,0.,0.), (100./255.,1.,1.),(1.,0.,0.)],
         'green' : [(0,0.,0.),(1.,0.,0.)] , 
         'blue' : [(0,0.,0.),(100./255.,0.,0.),(1.,1.,1.)] }
blue_red = LinearSegmentedColormap('BlueRed',cdict)

m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(-90,120,30), labels=[1,0,0,0])
m.drawmeridians(np.arange(-180,181,45),labels=[0,0,0,1])
m.imshow(np.flipud(snowcover),cmap=blue_red)
plt.title('MOD10C2: Eight Day Global Snow Cover')
plt.show()

此代码应显示积雪的图片。在

如果需要在不同的投影中处理数据，可以使用pythongdal接口将snowcover数组转换为地理位置数据集。在

将数据作为不规则网格处理也是可能的，但效率非常低。在

^{pr2}$

相应的经纬度网格。在