如何投影和重采样网格来匹配GDAL python的另一个网格？(How to project and

澄清：我不知何故离开了关键方面：不使用使用os.system或子 - 只要蟒蛇API。

我尝试转换NOAA GTX抵消电网的部分垂直基准转换，并没有完全跟随如何与蟒蛇为此在GDAL。我想花一格（在这种情况下，测深归因电网，但它可能是一个geotif），并用它作为我想要做的模板。如果我能做到这一点吧，我有一种感觉，这将极大地帮助人们利用这种类型的数据。

以下是我有一个肯定是行不通的。当我得到的目标数据集（dst_ds）运行gdalinfo，它不匹配源网格袋。

from osgeo import gdal, osr

bag = gdal.Open(bag_filename)
gtx = gdal.Open(gtx_filename)

bag_srs = osr.SpatialReference()
bag_srs.ImportFromWkt(bag.GetProjection())

vrt = gdal.AutoCreateWarpedVRT(gtx, None, bag_srs.ExportToWkt(), gdal.GRA_Bilinear,  0.125)

dst_ds = gdal.GetDriverByName('GTiff').Create(out_filename, bag.RasterXSize, bag.RasterYSize,
                                            1, gdalconst.GDT_Float32)
dst_ds.SetProjection(bag_srs.ExportToWkt())
dst_ds.SetGeoTransform(vrt.GetGeoTransform())

def warp_progress(pct, message, user_data):
  return 1

gdal.ReprojectImage(gtx, dst_ds, None, None, gdal.GRA_NearestNeighbour, 0, 0.125, warp_progress, None)

示例文件（但任何两个网格它们重叠，但在不同的投影会做）：

http://surveys.ngdc.noaa.gov/mgg/NOS/coast/F00001-F02000/F00574/BAG/ F00574_MB_2m_MLLW_2of3.bag
http://vdatum.noaa.gov/download/data/VDatum_National.zip MENHMAgome01_8301 / mllw.gtx

命令行相当于什么，我试图做的：

gdalwarp -tr 2 -2 -te 369179 4773093 372861 4775259 -of VRT -t_srs EPSG:2960 \
     MENHMAgome01_8301/mllw.gtx  mllw-2960-crop-resample.vrt
gdal_translate mllw-2960-crop-resample.{vrt,tif}

Answer 1:

由于杰米的答案。

#!/usr/bin/env python

from osgeo import gdal, gdalconst

# Source
src_filename = 'MENHMAgome01_8301/mllw.gtx'
src = gdal.Open(src_filename, gdalconst.GA_ReadOnly)
src_proj = src.GetProjection()
src_geotrans = src.GetGeoTransform()

# We want a section of source that matches this:
match_filename = 'F00574_MB_2m_MLLW_2of3.bag'
match_ds = gdal.Open(match_filename, gdalconst.GA_ReadOnly)
match_proj = match_ds.GetProjection()
match_geotrans = match_ds.GetGeoTransform()
wide = match_ds.RasterXSize
high = match_ds.RasterYSize

# Output / destination
dst_filename = 'F00574_MB_2m_MLLW_2of3_mllw_offset.tif'
dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Float32)
dst.SetGeoTransform( match_geotrans )
dst.SetProjection( match_proj)

# Do the work
gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)

del dst # Flush

Answer 2:

如果我理解正确的问题，你可以通过运行gdalwarp和gdal_translate作为子达成你的目标。只是组装你的选择，然后就比如下面：

import subprocess

param = ['gdalwarp',option1,option2...]
cmd = ' '.join(param)
process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout = ''.join(process.stdout.readlines())
stderr = ''.join(process.stderr.readlines())

if len(stderr) > 0:
    raise IOError(stderr)

它可能不是最完美的解决方案，但它会完成这项工作。一旦运行，只是用@中的数据加载到numpy的，并随身携带自己的方式。

文章来源: How to project and resample a grid to match another grid with GDAL python?