GDAL pythonで別のグリッドと一致するようにグリッドを投影して再サンプリングするにはどうすればいいですか？

Question

明確化：私は何とか重要な点を省いています：os.systemやサブプロセスを使用していません - 単にPython APIです。

私はNOAA GTXオフセットグリッドのセクションを垂直方向のデータ変換に変換しようとしていますが、これはGDALでpythonを使ってこれを行う方法に完全には従っていません。グリッド（この場合はBathimetry Attributed Gridですが、ジオテイフになる可能性があります）を使用して、テンプレートとして使用したいと思います。私がこのことを正しく行うことができれば、人々がこの種のデータを活用するのに大いに役立つと感じています。

これは私が間違いなく働いているものです。結果の宛先データセット（dst_ds）でgdalinfoを実行すると、ソースグリッドBAGと一致しません。

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) 

例ファイル（それらは重複の任意の2つのグリッドが、異なる投影しているだろう）：

• 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

答えのためのJamieに感謝します。

#!/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) 

これは最も洗練された解決策ではないかもしれませんが、それは完了します。それが実行されたら、gdalを使用してnumpyにデータをロードしてください。