MODIS AQUA数据-使用python GDAL堆叠/镶嵌数据
我知道如何使用gdal和python访问和绘制子数据集。然而，我想知道是否有一种方法可以使用包含在HDF4文件中的地理数据，以便我可以在许多年内查看同一区域。
如果可能的话，可以从数据中切出一个区域吗？
更新：
更具体地说:我绘制了MODIS数据，正如您在下面看到的那样，河流向下移动(矩形结构左上角)。因此，在一整年的时间里，它与我观察到的位置不同。
子数据集中有一个名为Geolocation Fields的目录，其中包含Long和Alt目录。那么，有没有可能访问这些信息或将其覆盖到数据上，以剪切出特定的区域？
例如，如果我们看一下美国宇航局的图片，是否有可能将其削减到10-15个大气压之间。和-5到0的长度。
您可以通过复制以下url来下载示例文件：
更新：
我跑了
x0, dx, dxdy, y0, dydx, dy = hdf_file.GetGeoTransform()这给了我以下输出：
x0: 0.0
dx: 1.0
dxdy: 0.0
y0: 0.0
dydx: 0.0
dy: 1.0以及
gdal.Warp(workdir2+"/output.tif",workdir1+"/MYD021KM.A2009002.1345.006.2012058153105.hdf")
这给了我以下错误：
ERROR 1: Input file /Volumes/Transcend/Master_Thesis/Data/AQUA_002_1345/MYD021KM.A2009002.1345.006.2012058153105.hdf has no raster bands.**更新2：**
下面是关于如何打开和读取hdf文件的代码：
all_files是一个包含如下文件名的列表：
MYD021KM.A2008002.1345.006.2012066153213.hdf
MYD021KM.A2008018.1345.006.2012066183305.hdf
MYD021KM.A2008034.1345.006.2012067035823.hdf
MYD021KM.A2008050.1345.006.2012067084421.hdf
etc .....
for fe in all_files:
print "\nopening file: ", fe
try:
hdf_file = gdal.Open(workdir1 + "/" + fe)
print "getting subdatasets..."
subDatasets = hdf_file.GetSubDatasets()
Emissiv_Bands = gdal.Open(subDatasets[2][0])
print "getting bands..."
Bands = Emissiv_Bands.ReadAsArray()
print "unit conversion ... "
get_name_tag = re.findall(".A(\d{7}).", all_files[i])[0]
print "name tag of current file: ", get_name_tag
# Code for 1 Band:
L_B_1 = radiance_scales[specific_band] * (Bands[specific_band] - radiance_offsets[specific_band]) # Source: MODIS Level 1B Product User's Guide Page 36 MOD_PR02 V6.1.12 (TERRA)/V6.1.15 (AQUA)
data_1_band['%s' % get_name_tag] = L_B_1
L_B_1_mean['%s' % get_name_tag] = L_B_1.mean()
# Code for many different Bands:
data_all_bands["%s" % get_name_tag] = []
for k in Band_nrs[lowest_band:highest_band]: # Bands 8-11
L_B = radiance_scales[k] * (Bands[k] - radiance_offsets[k]) # List with all bands
print "Appending mean value of {} for band {} out of {}".format(L_B.mean(), Band_nrs[k], len(Band_nrs))
data_all_bands['%s' % get_name_tag].append(L_B.mean()) # Mean radiance values
i=i+1
print "data added. Adding i+1 = ", i
except AttributeError:
print "\n*******************************"
print "Can't open file {}".format(workdir1 + "/" + fe)
print "Skipping this file..."
print "*******************************"
broken_files.append(workdir1 + "/" + fe)
i=i+1回答 1
Stack Overflow用户
发布于 2018-10-15 21:00:41
如果不知道确切的数据源和所需的输出等，就很难给你一个明确的答案。如上所述，您似乎拥有原生的.hdf格式的MODIS图像，并希望进行一些子集，以获得引用到相同区域的图像，然后绘制等。
从gdal模块查看gdal.Warp()可能会对您有所帮助。这种方法能够获取.hdf文件，并将一系列图像子集到具有相同分辨率/相同行数和列数的相同边界框中。
然后，您可以分析和绘制这些图像/比较像素等。
我希望这能给你一个良好的入门起点。
gdal.Warp文档：https://gdal.org/python/osgeo.gdal-module.html#Warp
更通用的扭曲帮助：https://www.gdal.org/gdalwarp.html
如下所示：
import gdal
# Set up the gdal.Warp options such as desired spatial resolution,
# resampling algorithm to use and output format.
# See: https://gdal.org/python/osgeo.gdal-module.html#WarpOptions
# for other options that can be specified.
warp_options = gdal.WarpOptions(format="GTiff",
outputBounds=[min_x, min_y, max_x, max_y],
xRes=res,
yRes=res,
# PROBABLY NEED TO SET t_srs TOO
)
# Apply the warp.
# (output_file, input_file, options)
gdal.Warp("/path/to/output_file.tif",
"/path/to/input_file.hdf",
options=warp_options)要编写的确切代码：
# Apply the warp.
# (output_file, input_file, options)
gdal.Warp('/path/to/output_file.tif',
'/path/to/HDF4_EOS:EOS_SWATH:"MYD021KM.A2009034.1345.006.2012058160107.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB',
options=warp_options)https://stackoverflow.com/questions/52763867
复制相似问题
腾讯云开发者
