Transforms a file with a vert datum crs to the wgs84 ellipsoid.

gdal_datum.py

from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import numpy as np
import blocks

ds = gdal.Open("n39_egm96.tif")
srs = ds.GetProjectionRef()

source = osr.SpatialReference()
source.ImportFromWkt(ds.GetProjection())

dest = osr.SpatialReference()
dest.ImportFromWkt(ds.GetProjection())
# Remove the vert CS so that it converts to the ellipsoid
dest.SetVertCS("")

transformer = osr.CoordinateTransformation(source, dest)

geotransform = ds.GetGeoTransform()

# Create the output file.
out_ds = gdal.GetDriverByName("GTiff").Create("n39_wgs84.tif", ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32)
gdal_array.CopyDatasetInfo(ds, out_ds)
out_band = out_ds.GetRasterBand(1)

band = ds.GetRasterBand(1)
for block in blocks.iter_blocks(band):
	print(block)
	data = np.array(blocks.read_block(band, block), dtype=float)
	for x in xrange(0, block[2]):
		c = x + block[0]
		for y in xrange(0, block[3]):
			r = y + block[1]
			geo_x,geo_y = gdal.ApplyGeoTransform(geotransform,c,r)			
			h = data[y,x]			
			dummyx, dummyy, new_h = transformer.TransformPoint(geo_x,geo_y,h)			
			data[y,x] = new_h
			blocks.write_block(out_band, block, data)