gdal_grid with CSV file and VRT "No point geometry found on layer"

Question

Given the following .vrt file which I manually created:

<OGRVRTDataSource>
    <OGRVRTLayer name="DGM5_BE">
        <SrcDataSource>DGM5_BE.txt</SrcDataSource>
        <SrcLayer>DGM5_BE</SrcLayer>
        <LayerSRS>EPSG:25833</LayerSRS>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="x" y="y" z="d"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

and the associated DGM5_BE.txt - here are the first lines:

x;y;d
397200;5837250;59,15
397205;5837250;59,04
397200;5837245;59,03
397205;5837245;58,94
397195;5837240;58,93
397200;5837240;58,9
397205;5837240;58,86
397210;5837240;58,9
397190;5837235;58,98

I run the following command to convert the DGM5:

$ gdal_grid -of GTiff -ot Float32 -l DGM5_BE DGM5_BE.vrt DGM5_BE.tif

This fails with the following error message:

ERROR 1: Failed to open datasource `DGM5_BE.txt'.
No point geometry found on layer DGM5_BE, skipping.

Answer 1

From the CSV driver documentation:

Starting with GDAL 1.8.0, for files structured as CSV, but not ending with .CSV extension, the 'CSV:' prefix can be added before the filename to force loading by the CSV driver.

Either rename DGM5_BE.txt to DGM5_BE.csv or change the <SrcDataSource> element to:

<SrcDataSource>CSV:DGM5_BE.txt</SrcDataSource>

Answer 2

For what its worth I created a demo that shows example usage of gdal_grid.

def main():
    """
    Demo using gdal_grid

    References:
        https://gis.stackexchange.com/questions/254330/python-gdal-grid-correct-use
        https://gdal.org/programs/gdal_grid.html

    """
    from os.path import exists
    import ubelt as ub
    import geopandas as gpd
    import numpy as np

    # Generate Demo Data Ungridded data
    # Generate a random value at points along the African continent
    xyz = []
    world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
    rng = np.random.RandomState(0)
    for shape in world[world['continent'] == 'Africa'].geometry:
        if shape.type == 'MultiPolygon':
            for geom in shape.geoms:
                for x, y in geom.exterior.coords:
                    z = rng.rand()
                    xyz.append((x, y, z))
        elif shape.type == 'Polygon':
            for x, y in shape.exterior.coords:
                z = rng.rand()
                xyz.append((x, y, z))
    xyz = np.array(xyz)

    # Write the demo data to a CSV file
    lines = ['x;y;d']
    for x, y, z in xyz:
        lines.append(f'{x};{y};{z}')
    text = '\n'.join(lines)
    with open('tmp.csv', 'w') as file:
        file.write(text)

    # I was unable to get the command to work with the CSV file alone
    # I needed to create a VRT that pointed at the CVS file.
    with open('tmp.vrt', 'w') as file:
        file.write(ub.codeblock(
            '''
            <OGRVRTDataSource>
                <OGRVRTLayer name="tmp">
                    <SrcDataSource>tmp.csv</SrcDataSource>
                    <SrcLayer>tmp</SrcLayer>
                    <LayerSRS>EPSG:4326</LayerSRS>
                    <GeometryType>wkbPoint</GeometryType>
                    <GeometryField encoding="PointFromColumns" x="x" y="y" z="d"/>
                </OGRVRTLayer>
            </OGRVRTDataSource>
            '''))

    # Can modify these but they are not needed
    minx, miny = xyz.T[0:2].min(axis=1)
    maxx, maxy = xyz.T[0:2].max(axis=1)
    lonext = maxx - minx
    latext = maxy - miny
    ar = latext / lonext
    xres = int(512)
    yres = int(xres * ar)

    # Call the GDAL Grid command
    command = (
        'gdal_grid '
        # '-ot Float32 -of GTiff '
        # '-zfield d '
        # '-a_srs EPSG:4326 '
        # f'-txe {minx - 1} {maxy + 1} '
        # f'-tye {miny - 1} {maxx + 1} '
        # f'-outsize {xres} {yres} '
        # '-a invdist:power=2.0:smoothing=1.0 '
        # '-a nearest:radius1=1.0:radius2=1.0 '
        # '-a linear:radius=1.0 '
        'tmp.vrt tmp.tif')

    ub.delete('tmp.tif')
    print(command)
    info = ub.cmd(command, verbose=3)
    assert info['ret'] == 0

    if exists('tmp.tif'):
        # Show the result if it worked
        import tifffile
        data = tifffile.imread('tmp.tif')
        from matplotlib import pyplot as plt
        print('data.shape = {!r}'.format(data.shape))
        plt.imshow(data)

The result is this: