Search code examples
pythonpyvista

How to read dem files in python pyvista


i have .dem files like those here: http://ddfe.curtin.edu.au/models/ERTM2160/data/dem/

In python pyvista i have eg:

import pyvista as pv
file = 'pick_one_from_the_link_above.dem'
mesh = pv.read(file)

The output says:

mesh.dimensions
[-2147483648,-2147483648,1]

which ist, except for the minus sign, the square root of mesh.n_points

Trying to plot or extract the points with mesh.points i get a message that negative dimensions are not allowed. Trying the folowing:

mesh.dimensions = [int(numpy.sqrt(mesh.n_points)),int(numpy.sqrt(mesh.n_points)),1]

results in an error message:

OverflowError: SetDimensions argument 1: value is out of range for int

Can someone tell me what i am doing wrong, i have no idea? Or maybe has an idea how to read those files for making a surface plot?

Thanks alot :)


Solution

  • @larsks is correct in the comment above. These ".dem" files are not in the format that PyVista and the VTK reader it wraps expect. You should use np.fromfile to read the data: arr = np.fromfile('N00E015.dem', dtype=np.int16). Further from the docs listed at your link:

    Due to its total size of 44 GB, the model is partitioned and distributed in terms of 881 binary files of 5 deg x 5 deg size for each functional. Each 5 deg x 5 deg tile contains 2500 x 2500 grid points in centre-of-cell representation (grid points are not situated on integer meridians and parallels)

    You'll just create a pv.UniformGrid of that size and add the data. For instance:

    import numpy as np
    import pyvista as pv
    
    arr = np.fromfile('N00E015.dem', dtype=np.int16)
    
    grid = pv.UniformGrid()
    grid.dimensions = (2500, 2500, 1)
    grid.origin = (0, 0, 0) # you need to figure this out
    grid['dem'] = arr
    
    grid.plot()
    

    enter image description here

    In order to get the correct spatial reference for the grid, you'll need to set the origin point of each subset/grid.

    Further, the PyVista community is more active on the PyVista support forum than on Stack Overflow: https://github.com/pyvista/pyvista-support