search

kylebarron / demquery

Wrapper around rasterio to query points on a Digital Elevation Model
MIT License
11 stars 0 forks source link

AssertionError when trying to use Netherlands DEMs #3

Open craigxchen opened 3 years ago

craigxchen commented 3 years ago

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
kylebarron commented 3 years ago

Can you create an minimal example of using demquery directly? E.g. try using query_points(points) directly and find the point where it fails?

craigxchen commented 3 years ago

Hi Kyle,

Thanks for responding.

I tried it using a minimal query and it's failing on every single point of the graph.

It seems like the assertion error doesn't have anything to do with the specific points right? It's coming from the _check_bounds method. My best guess is that somehow the TIFs aren't getting tiled together correctly?

kylebarron commented 3 years ago

I think it's likely your DEM is not in the same projection as your input points. This package won't do any automatic reprojection of the raster file, and if the raster data is in some local projection, it won't work and _check_bounds will fail.

craigxchen commented 3 years ago

Do you have any suggestions on how I can fix the projection issue? Do you think there's a quick fix or is it something more involved?

kylebarron commented 3 years ago

You can look into GDAL/rasterio for reprojecting your data. https://gdal.org/programs/gdalwarp.html