3

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.

3

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>
2
  • Thx. I can't believe how long it takes to calculate a TIFF! I already added --config GDAL_NUM_THREADS ALL_CPUS. Is there recommendable format other than TIFF which can be used in QGIS?
    – JJD
    Jul 4 '14 at 13:47
  • The problem is that triangulation is only a single thread process, irrespective of how many CPU threads you allow this process it can only use 1 to triangulate with the other threads for reading/writing. There are options that do slow down a TIFF, what compression are you using? I don't think you will get any more performance from this process regardless of the output format. Aug 9 '18 at 3:58
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:

enter image description here

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.