kylebarron / demquery

Wrapper around rasterio to query points on a Digital Elevation Model

demquery

Wrapper around rasterio to query points on a Digital Elevation Model.

Features

• Use multiple raster files without having to merge them into a new file
• Query many points at once
• Optional 2D interpolation (linear, cubic, or quintic)
• Reasonably performant by reading the minimal data required from raster

Install

I recommend first installing dependencies with Conda, then installing demquery itself with pip.

conda install gdal rasterio numpy scipy -c conda-forge

pip install demquery

CLI Script

> demquery --help
Usage: demquery [OPTIONS] FEATURES...

  Assign elevations to GeoJSON

Options:
  -d, --dem PATH          Paths to DEM files.  [required]
  -g, --dem-glob TEXT     Glob expression for DEM paths if folder is provided.
  -b, --band INTEGER      Band of rasters to use  [default: 1]
  -i, --interp-kind TEXT  either None, "linear", "cubic", or "quintic". None
                          will do no interpolation and choose the value in the
                          DEM closest to the provided point. linear creates a
                          3x3 grid and runs linear interpolation; cubic
                          creates a 5x5 grid and runs cubic interpolation;
                          quintic creates a 7x7 grid and runs quintic
                          interpolation.
  --help                  Show this message and exit.

echo \
    '{"type":"Feature","properties":{"name": "Glacier Peak"},"geometry":{"type":"Point","coordinates":[-121.2436843,48.0163834]}}' \
    | demquery -d /path/to/dem/files

Outputs:

{"type": "FeatureCollection", "features": [{"type": "Feature", "geometry": {"type": "Point", "coordinates": [-121.243684, 48.016383, 1431.5755615234375]}, "properties": {"name": "Glacier Peak"}}]}

Documentation

from demquery import Query

dem_paths = ['dem1.tif', 'dem2.tif']
query = Query(dem_paths)

# Points must be in longitude, latitude order!
# These points are in Manhattan, not Antarctica
points = [(-73.985564, 40.757965), (-73.968520, 40.778912)]
elevations = query.query_points(points, interp_kind='linear')

Data Download

For a great visual tool to download worldwide SRTM data, check out these sites:

• 30m resolution: http://dwtkns.com/srtm30m/
• 90m resolution: http://dwtkns.com/srtm/

Releasing

To upload a new release to PyPI

python setup.py sdist
twine upload dist/demquery-0.3.0.tar.gz

demquery's Issues

Update conda version

Conda is still on 0.2.0 I believe

AssertionError when trying to use Netherlands DEMs

• demquery version: 0.3.1
• Python version: 3.8.6
• Operating System: Windows 10

Description

Trying to add elevation values to a map of Amsterdam using TIFs from pdok.nl. I'm confident that the TIF files actually cover all of the graph.

What I Did

Here is the code I ran:

dem_paths = glob.glob("../data/netherlands_dem_5m/*.tif")

graph = ox.graph.graph_from_place(name, network_type="walk")
test_elev = add_node_elevations_from_dem(graph.nodes, dem_paths, interp_kind="linear")
t = dict((k, test_elev[k]['elevation']) for k in test_elev.keys())
nx.set_node_attributes(graph, name="elevation", values=t)
ox.save_graphml(G=graph, filepath=f"../data/{name}.graphml")

And the error message:

  File "C:\Users\craig\Documents\Embarc\Coding\elevation_injector\src\netherlands_query.py", line 22, in <module>
    test_elev = add_node_elevations_from_dem(graph.nodes, dem_paths, interp_kind="linear")

  File "C:\Users\craig\Documents\Embarc\Coding\elevation_injector\src\dem_elevation.py", line 32, in add_node_elevations_from_dem
    elevations = query.query_points(points)

  File "C:\Users\craig\Anaconda3\envs\embarc\lib\site-packages\demquery\demquery.py", line 74, in query_points
    self._check_bounds(dem, points, num_buffer=num_buffer)

  File "C:\Users\craig\Anaconda3\envs\embarc\lib\site-packages\demquery\demquery.py", line 150, in _check_bounds
    assert maxrow <= dem.height

AssertionError

Add more tests

Add a couple tests with low resolution DEM data. Not sure what the lowest-resolution dem would be. An SRTM tile from here is 35MB zipped.

Something that size is probably a good size to download (and cache) for tests, but not to include in the repository.