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

Aclaración: De alguna manera, dejé de lado 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 desplazamiento NOAA GTX para transformaciones de referencia vertical y no estoy siguiendo totalmente cómo hacerlo en GDAL con python. Me gustaría tomar una cuadrícula (en este caso, una Cuadrícula Atribuida de Batimetría, pero podría ser un geotif) y usarla como la plantilla a la que me gustaría hacer. Si puedo hacer esto bien, tengo la sensación de que ayudará mucho a las personas a utilizar 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) 

Los archivos de ejemplo (pero cualquiera de las dos cuadrículas donde se superponen, pero están en diferentes proyecciones harían):

  • 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 equivalente 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} 

    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 

    Si entiendo la pregunta correctamente, podría lograr su objective ejecutando gdalwarp y gdal_translate como subprocesos. Simplemente reúna sus opciones y luego 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 que no sea la solución más elegante, pero hará el trabajo. Una vez que se ejecuta, simplemente cargue sus datos en numpy usando gdal y siga su camino.