我目前在R
中有一个处理链,它下载MODIS数据,然后从系统调用gdalwarp
重新投影特定的子数据集(例如NDVI)到WGS1984。然后将生成的GeoTiffs
收集到HDF5
文件中进行进一步处理。
现在我正在将处理链移植到python
,想知道是否有一种方法可以使用gdal
模块的功能跳过将GeoTiffs
写入磁盘的步骤。
明确一下问题:
我是否能够仅使用gdal
模块的python绑定执行gdalwarp
而不写入磁盘?
我进行了一些研究,最接近我的问题的答案是这些帖子:
第一种方法需要一个模板,所以不是我要找的。
第二种方法看起来更有前途,它使用了AutoCreateWarpedVRT
函数,这似乎是我想要的。虽然与答案中的示例相反,我的结果与参考结果不匹配(与任何误差阈值无关)。
在直接调用gdalwarp
的以前实现中,除了目标参考系统外,我还指定了目标分辨率。因此,我认为这可能是有所区别的 - 但我无法在python中通过gdal绑定设置它。
以下是我的尝试(抱歉,没有MODIS数据无法重现):
import gdal
import osr
ds = gdal.Open('/data/MOD13A2.A2016305.h18v07.005.2016322013359.hdf')
t_srs = osr.SpatialReference()
t_srs.ImportFromEPSG(4326)
src_ds = gdal.Open(ds.GetSubDatasets()[0][0], gdal.GA_ReadOnly)
dst_wkt =t_srs.ExportToWkt()
error_threshold = 0.125
resampling=gdal.GRA_NearestNeighbour
tmp_ds = gdal.AutoCreateWarpedVRT( src_ds,
None, # src_wkt : left to default value --> will use the one from source
dst_wkt,
resampling,
error_threshold)
# create tiff
dst_ds = gdal.GetDriverByName('GTiff').CreateCopy('warp_test.tif', tmp_ds)
dst_ds = None
这是供参考的:
gdalwarp -ot Int16 -tr 0.00892857142857143 0.00892857142857143 -t_srs EPSG:4326 "HDF4_EOS:EOS_GRID:MOD13A2.A2016305.h18v07.005.2016322013359.hdf:MODIS_Grid_16DAY_1km_VI:1 km 16 days NDVI" MOD13A2.A2016305.h18v07.005.2016322013359_MODIS_Grid_16DAY_1km_VI_1_km_16_days_NDVI.tif
比较:
i1 = gdal.Open('warp_test.tif')
i2 = gdal.Open('MOD13A2.A2016305.h18v07.005.2016322013359_MODIS_Grid_16DAY_1km_VI_1_km_16_days_NDVI.tif')
# test
print(i1.RasterXSize,i1.RasterYSize)
1267 1191
#reference
print(i2.RasterXSize,i2.RasterYSize)
1192 1120
i1.GetRasterBand(1).Checksum() == i2.GetRasterBand(1).Checksum()
False
所以您可以看到,使用gdal.AutoCreateWarpedVRT
函数会导致数据集具有不同的尺寸和分辨率。
sr
、dst_wkt
和tmp_ds2
是什么?它们无法被重现,因此很难理解问题是什么。 - Rutger Kassies