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}
1
  • 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 8, 2012 at 16:24

2 Answers 2

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
4
  • Works perfectly - you could wrap it up into a script that gets the filenames from sys.argv for neatness.
    – Spacedman
    May 11, 2012 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, 2020 at 19:39
  • 1
    this works great, just a detail: local variable 'src_geotrans' is assigned to but never used
    – duff18
    Feb 16, 2022 at 7: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, 2022 at 15:25
4

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.

3
  • subprocess would definitely get the image corrected. However, I somehow missed mentioning that I was aiming for a GDAL python API only solution. May 10, 2012 at 16:16
  • 2
    Would be nice to see the options required to do what the python script does anyway.
    – Spacedman
    May 11, 2012 at 12:28
  • The library version of the app programs are available via the swig based python wrappers. May 1, 2020 at 13:57

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.