1

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?

  • I recommend to ask this from the gdal-dev mailing list. There are some other posts about vertical datums also lists.osgeo.org/pipermail/gdal-dev/2019-May/date.html. I also recommend to update to GDAL 3 because of big changes with projections trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn. – user30184 Jun 5 '19 at 8:12
  • 1
    I don't know if this has any bearing on your problem, but the UTM coordinates are probably swapped. The 9 million value is likely the northing / Y while the 787k is the easting/X. – mkennedy Jun 5 '19 at 20:07
  • Please make an answer from the feedback that you got from the mailing list. – user30184 Jun 6 '19 at 8:06
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

| improve this answer | |
  • Please continue on the mailing list if GDAL 3 did not give the right answer for you. – user30184 Jun 10 '19 at 14:37
1

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.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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