¿Cómo proyectar y volver a muestrear una cuadrícula para que coincida con otra cuadrícula con GDAL python?

Question

Aclaración: de alguna manera omití el aspecto clave: no usar os.system o subprocess, solo la API de python.

Estoy tratando de convertir una sección de una cuadrícula de compensación NOAA GTX para transformaciones verticales de datum y no seguir totalmente cómo hacer esto en GDAL con python. Me gustaría tomar una grilla (en este caso, una Grilla Atribuida de Batimetría, pero podría ser un geotif) y usarla como la plantilla que me gustaría hacer. Si puedo hacer esto bien, tengo la sensación de que será de gran ayuda que las personas hagan uso de este tipo de datos.

Esto es lo que tengo que definitivamente no funciona. Cuando ejecuto gdalinfo en el conjunto de datos de destino resultante (dst_ds), no coincide con la BAG de la cuadrícula de origen.

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) 

archivos de ejemplo (pero ningún dos rejillas donde se superponen, pero se encuentran en diferentes proyecciones haría):

• 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

La línea de comando equivale a lo que estoy tratando de hacer:

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

Gracias a Jamie por la respuesta.

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

Si entiendo la pregunta correctamente, puede lograr su objetivo ejecutando gdalwarp y gdal_translate como subprocesos. Sólo montar sus opciones a continuación, haga lo siguiente, por ejemplo:

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) 

puede no ser la solución más elegante, pero va a hacer el trabajo. Una vez que se ejecuta, simplemente cargue sus datos en numpy usando gdal y siga su camino.