17

Clarification: I somehow left out the key aspect: not using os.system or subprocess - just the python API.

I'm trying to convert a section of a NOAA GTX offset grid for vertical datum transformations and not totally following how to do this in GDAL with python. I'd like to take a grid (in this case a Bathymetry Attributed Grid, but it could be a geotif) and use it as the template that I'd like to do to. If I can do this right, I have a feeling that it will greatly help people make use of this type of data.

Here is what I have that is definitely not working. When I run gdalinfo on the resulting destination dataset (dst_ds), it does not match the source grid 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)

Example files (but any two grids where they overlap, but are in different projections would do):

The command line equivalent to what I'm trying to do:

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}
Wraf
  • 747
  • 2
  • 10
  • 24
Kurt Schwehr
  • 2,638
  • 3
  • 24
  • 41
  • What is the output of the WKT for bag_srs? Have you verified it is the same SRS that "bag" gives? I've come across some WKT that is... well, not well formatted... I notice the command line version specifies EPSG:2960 (which is NAD83?). I haven't used gdal in a long while, but if I ran into this I guess I would start by making sure that the reprojection is using the proper SRS values. Sorry, not a good answer ... that's why it's a comment. – Kasapo May 08 '12 at 16:24

2 Answers2

31

Thanks to Jamie for the answer.

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

# 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
Kurt Schwehr
  • 2,638
  • 3
  • 24
  • 41
  • Works perfectly - you could wrap it up into a script that gets the filenames from sys.argv for neatness. – Spacedman May 11 '12 at 12:27
  • what if the source data and match data only have part overlap? Any idea about set the nodata value for non-overlapping areas? – Xue Apr 29 '20 at 19:39
  • 1
    this works great, just a detail: local variable 'src_geotrans' is assigned to but never used – duff18 Feb 16 '22 at 07:52
  • Is there a command line gdal equivalent to this? similar to the original post but reading `-te` from the match_file instead of having to manually specify – konstanze Mar 25 '22 at 15:25
3

If I understand the question correctly, you could accomplish your goal by running gdalwarp and gdal_translate as subprocesses. Just assemble your options then do the following for example:

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)

It may not be the most elegant solution, but it will get the job done. Once it is run, just load your data into numpy using gdal and carry on your way.

Nathan Smith
  • 118
  • 1
  • 6