我目前有一个加工链R
下载MODIS数据然后调用gdalwarp
从系统将特定子数据集(例如 NDVI)重新投影到 WGS1984 中。所结果的GeoTiffs
然后被收集到一个HDF5
文件以供进一步处理。
现在我将处理链移至python
,我想知道是否有办法跳过写作步骤GeoTiffs
到磁盘的功能gdal
module.
需要明确的是,问题是:
我可以表演吗gdalwarp
严格使用来自的 python 绑定gdal
模块并且不写入磁盘?
我做了一些研究,最接近我的问题的答案是这些帖子:
- 如何使用 GDAL python 投影和重新采样网格以匹配另一个网格? https://stackoverflow.com/questions/10454316/how-to-project-and-resample-a-grid-to-match-another-grid-with-gdal-python
- 使用 gdal Python 绑定复制 gdalwarp 的结果 https://gis.stackexchange.com/questions/139906/replicating-result-of-gdalwarp-using-gdal-python-bindings
第一种方法需要模板,所以不是我真正想要的。
第二种方法看起来更有希望,它使用该函数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
函数会产生具有不同维度和分辨率的数据集。