来自 hdf 文件的经纬度信息 python
lat,lon information from hdf file python
我有一个 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'))[:]
它给我一个错误：
HDF4Error: select: non-existent dataset
我试过：
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
您使用的数据文件是MODIS Level 3 产品。所有级别 3 的产品都被插值到一些规则的网格上。在 MOD10C2 的情况下，网格是所谓的气候建模网格 (CMG)。该网格规则地间隔 0.05 度。 Panoply 知道这一点。
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()
此代码应显示积雪的图片。
如果您需要使用不同投影中的数据，您可以使用 python GDAL 接口将 snowcover 数组转换为地理定位数据集。
也可以将数据作为不规则网格处理，但效率很低。
lon,lat = np.meshgrid(np.arange(-180,180,0.05),np.arange(-90,90,0.05))
将是相应的经度和纬度网格。
通常情况下，纬度和经度信息不在 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()
如果您不知道确切的变量名称，您可以尝试 pyhdf.V 使用以下显示包含在其中的 vgroup 内容的程序
任何 HDF 文件。
from pyhdf.HDF import *
from pyhdf.V import *
from pyhdf.VS import *
from pyhdf.SD import *
def describevg(refnum):
# Describe the vgroup with the given refnum.
# Open vgroup in read mode.
vg = v.attach(refnum)
print "----------------"
print "name:", vg._name, "class:",vg._class, "tag,ref:",
print vg._tag, vg._refnum
# Show the number of members of each main object type.
print "members: ", vg._nmembers,
print "datasets:", vg.nrefs(HC.DFTAG_NDG),
print "vdatas: ", vg.nrefs(HC.DFTAG_VH),
print "vgroups: ", vg.nrefs(HC.DFTAG_VG)
# Read the contents of the vgroup.
members = vg.tagrefs()
# Display info about each member.
index = -1
for tag, ref in members:
index += 1
print "member index", index
# Vdata tag
if tag == HC.DFTAG_VH:
vd = vs.attach(ref)
nrecs, intmode, fields, size, name = vd.inquire()
print " vdata:",name, "tag,ref:",tag, ref
print " fields:",fields
print " nrecs:",nrecs
vd.detach()
# SDS tag
elif tag == HC.DFTAG_NDG:
sds = sd.select(sd.reftoindex(ref))
name, rank, dims, type, nattrs = sds.info()
print " dataset:",name, "tag,ref:", tag, ref
print " dims:",dims
print " type:",type
sds.endaccess()
# VS tag
elif tag == HC.DFTAG_VG:
vg0 = v.attach(ref)
print " vgroup:", vg0._name, "tag,ref:", tag, ref
vg0.detach()
# Unhandled tag
else:
print "unhandled tag,ref",tag,ref
# Close vgroup
vg.detach()
# Open HDF file in readonly mode.
filename = 'yourfile.hdf'
hdf = HDF(filename)
# Initialize the SD, V and VS interfaces on the file.
sd = SD(filename)
vs = hdf.vstart()
v = hdf.vgstart()
# Scan all vgroups in the file.
ref = -1
while 1:
try:
ref = v.getid(ref)
print ref
except HDF4Error,msg: # no more vgroup
break
describevg(ref)
我认为问题在于该文件没有传统的纬度和经度数据（与许多 .nc 文件一样）。
当我想处理 MYD14 数据时遇到类似的问题（这是一个关于 fire-mask 的 MODIS 文件）。我搜索了很长时间才能找到解决方案。这是我的发现：
①如果MODIS文件使用SIN-Grid(Sinusoidal Projection)定义数据，文件不会给你传统的lon,lat数据。
②更多Sinusoidal Projection的细节，可以看这个网址：https://code.env.duke.edu/projects/mget/wiki/SinusoidalMODIS
我有一个 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'))[:]
它给我一个错误：
HDF4Error: select: non-existent dataset
我试过：
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
您使用的数据文件是MODIS Level 3 产品。所有级别 3 的产品都被插值到一些规则的网格上。在 MOD10C2 的情况下，网格是所谓的气候建模网格 (CMG)。该网格规则地间隔 0.05 度。 Panoply 知道这一点。
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()
此代码应显示积雪的图片。
如果您需要使用不同投影中的数据，您可以使用 python GDAL 接口将 snowcover 数组转换为地理定位数据集。
也可以将数据作为不规则网格处理，但效率很低。
lon,lat = np.meshgrid(np.arange(-180,180,0.05),np.arange(-90,90,0.05))
将是相应的经度和纬度网格。
通常情况下，纬度和经度信息不在 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()
如果您不知道确切的变量名称，您可以尝试 pyhdf.V 使用以下显示包含在其中的 vgroup 内容的程序 任何 HDF 文件。
from pyhdf.HDF import *
from pyhdf.V import *
from pyhdf.VS import *
from pyhdf.SD import *
def describevg(refnum):
# Describe the vgroup with the given refnum.
# Open vgroup in read mode.
vg = v.attach(refnum)
print "----------------"
print "name:", vg._name, "class:",vg._class, "tag,ref:",
print vg._tag, vg._refnum
# Show the number of members of each main object type.
print "members: ", vg._nmembers,
print "datasets:", vg.nrefs(HC.DFTAG_NDG),
print "vdatas: ", vg.nrefs(HC.DFTAG_VH),
print "vgroups: ", vg.nrefs(HC.DFTAG_VG)
# Read the contents of the vgroup.
members = vg.tagrefs()
# Display info about each member.
index = -1
for tag, ref in members:
index += 1
print "member index", index
# Vdata tag
if tag == HC.DFTAG_VH:
vd = vs.attach(ref)
nrecs, intmode, fields, size, name = vd.inquire()
print " vdata:",name, "tag,ref:",tag, ref
print " fields:",fields
print " nrecs:",nrecs
vd.detach()
# SDS tag
elif tag == HC.DFTAG_NDG:
sds = sd.select(sd.reftoindex(ref))
name, rank, dims, type, nattrs = sds.info()
print " dataset:",name, "tag,ref:", tag, ref
print " dims:",dims
print " type:",type
sds.endaccess()
# VS tag
elif tag == HC.DFTAG_VG:
vg0 = v.attach(ref)
print " vgroup:", vg0._name, "tag,ref:", tag, ref
vg0.detach()
# Unhandled tag
else:
print "unhandled tag,ref",tag,ref
# Close vgroup
vg.detach()
# Open HDF file in readonly mode.
filename = 'yourfile.hdf'
hdf = HDF(filename)
# Initialize the SD, V and VS interfaces on the file.
sd = SD(filename)
vs = hdf.vstart()
v = hdf.vgstart()
# Scan all vgroups in the file.
ref = -1
while 1:
try:
ref = v.getid(ref)
print ref
except HDF4Error,msg: # no more vgroup
break
describevg(ref)
我认为问题在于该文件没有传统的纬度和经度数据（与许多 .nc 文件一样）。 当我想处理 MYD14 数据时遇到类似的问题（这是一个关于 fire-mask 的 MODIS 文件）。我搜索了很长时间才能找到解决方案。这是我的发现： ①如果MODIS文件使用SIN-Grid(Sinusoidal Projection)定义数据，文件不会给你传统的lon,lat数据。 ②更多Sinusoidal Projection的细节，可以看这个网址：https://code.env.duke.edu/projects/mget/wiki/SinusoidalMODIS