How to Set the Vertical Coordinate system when creating 3D Shapefile?

Question

I have data of wells in CSV format, which I'm converting to a 3D shapefile, using ogr2ogr. The latitude and longitude are in WGS84, while the depth is from the WGS84 geoid. I have the OSGeo4W installation of ogr2ogr.

The process completes successfully, and Points with XYZ are created;

However I can't figure out how to set the Vertical Coordinate system when creating this shapefile.

My CSV file is like this:

id, name, Lat, Long, depth
1, A, 12, 35, -10
2, B, 13, 40, -15
3, C, 15, 45, -20

My vrt file is like this:

<OGRVRTDataSource>
    <OGRVRTLayer name="data">
        <SrcDataSource>data.csv</SrcDataSource>
        <GeometryType>wkbPoint25D</GeometryType>
        <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="Long" y="Lat" z="depth"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

I'm using the following command to create the shapefile:

ogr2ogr -f "ESRI Shapefile" input.shp data.vrt

How to Set the Vertical Coordinate system when creating 3D Shapefile with ogr2ogr?

Answer 1

I believe you can specify a compound coordinate system (vertical + horizontal) as "EPSG:X+Y" where X and Y are the EPSG codes for your horizontal and vertical coordinate systems.

The following reports correctly in ogrinfo and converts with ogr2ogr:

<OGRVRTDataSource>
    <OGRVRTLayer name="data">
        <SrcDataSource>data.csv</SrcDataSource>
        <GeometryType>wkbPoint25D</GeometryType>
        <LayerSRS>EPSG:4326+5703</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="Long" y="Lat" z="depth"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

This entry in the gdal mailing list is helpful.

For whatever reason, that doesn't add the vertical coordinate system into the prj file when using ogr2ogr. In order to get the prj output correctly, you could either update/append into an existing (empty) shapefile, or force the .prj wkt from another wkt string, such as:

ogr2ogr -t_srs 'path/to/myWKT.wkt' output.shp myVRT.vrt

You can change it to 'ESRI::myWKT.wkt' if you want to push an ESRI WKT.