Transform datum CSV->CSV using ogr2ogr

Question

First I have a set of points in UTM 48S / WGS84 (EPSG 32748) as test.csv:

X,Y,Z
9231824.04,787648.04,711.1442

Second I defined a test.vrt for ogr2ogr to read the data as follows

<OGRVRTDataSource>
    <OGRVRTLayer name="test">
        <SrcDataSource>test.csv</SrcDataSource>
        <GeometryType>wkbPoint25D</GeometryType>
        <LayerSRS>EPSG:32748+4979</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

Thirdly I then convert to UTM 48S using EGM2008 datum as follows:

ogr2ogr -f CSV output.csv test.vrt -t_srs EPSG:32748+3855

The resulting file output.csv equals though the source:

X,Y,Z
9231824.04,787648.04,711.1442

I would expect for the Z-value to be difference. My OGR version is as follow:

ogr2ogr --version
GDAL 2.2.3, released 2017/11/20

What am I doing wrong?

Answer 1

On the GDAL-DEV mailing list I was kindly adressed to install vertical datum grids for GDAL. I did not achieve this under Windows using installers from gisinternal. On Ubuntu though I followed hints from https://gis.stackexchange.com/a/258769/144133 and copied http://download.osgeo.org/proj/vdatum/egm08_25/egm08_25.gtx to /usr/share/proj/. The installation was tested as follows

echo "9231824.04 787648.04 711.1442" | gdaltransform -s_srs EPSG:32748  -t_srs EPSG:32748+3855

giving the answer

9231824.03999554 787648.040034059 689.737249081086

For GDAL v3 I was recommended to first convert in 2D to WGS84 and then convert datum to EGM2008 and reproject. The following is untested.

echo "9231824.04 787648.04 711.1442" | \
   gdaltransform -s_srs EPSG:32748 -t_srs EPSG:4326 | \
   gdaltransform -s_srs EPSG:4979  -t_srs EPSG:32748+3855

When running now

ogr2ogr -f CSV output.csv test.vrt -s_srs EPSG:32748+4326 -t_srs EPSG:32748+3855

The answer is still unchanged 711.1442

Answer 2

Since I could not get ogr2ogr to work using a VRT-file I used Python to loop through a CSV file using gdaltransform. The procedure is as follows:

Source CSV format for data_wgs84.csv:

id,E,N,Z,Lat,Lon
G2259MH005,784557.8642,9233889.868,729.9722,-6.923923467,107.5750634
G2259MH005,784560.1291,9233903.801,729.964,-6.923797448,107.5750832
G2606MH001,787774.9186,9235390.141,737.2705,-6.910207604,107.604081
G2606MH002,787791.0091,9235467.886,738.0887,-6.909504258,107.6042226

The Python snipped is using Pandas library as follows

import subprocess
import pandas as pd 

file = r'data_wgs84.csv'
df = pd.read_csv(file)

for index, row in df.iterrows():
    cmd = 'echo {} {} {} | gdaltransform -s_srs EPSG:32748 -t_srs EPSG:32748+3855'.format(row['E'], row['N'], row['Z'])
    output = subprocess.check_output(cmd, shell=True)
    df.at[index,'E'] = (float(output.split()[0]))
    df.at[index,'N'] = (float(output.split()[1]))
    df.at[index,'Z'] = (float(output.split()[2]))

df.to_csv('data_egm2008.csv', encoding='utf-8', index=False

The output format of the data_egm2008.csv is reading:

id,E,N,Z,Lat,Lon
G2259MH005,784557.864199999,9233889.868,708.093566678607,-6.923923467000001,107.5750634
G2259MH005,784560.129100001,9233903.801,708.085797834495,-6.923797447999999,107.5750832
G2606MH001,787774.9186,9235390.141,715.342885635243,-6.910207604,107.60408100000001
G2606MH002,787791.0091,9235467.886,716.161568373874,-6.909504257999999,107.6042226

I suppose there's room for improvement.