Recommended workflow for map combining topography and street network

Specific Domains
,
1

I want to plot what I imagine (perhaps naively) to be a simple map, combining a topographic layer and the street network, for a small geographic area (one city district). I believe I cannot use an already-existing layer combining streets and topography: I want to magnify the small variations in elevation within a city district. I have managed to do a subset of this task (topographic map) with GMT.jl but I find the documentation difficult to navigate, and so it is mostly trial and error attempting to combine examples that have nothing to do with each other and use different subsets of GMT. Can anyone recommend a tool or workflow for this task?

1 Like
2

Hi @fa-bien , try to share the data you want to plot. People can then share examples with different packages.

Below is a random example with GeoStats.jl + GeoArtifacts.jl:

using GeoStats
using GeoArtifacts

import GLMakie as Mke

# download artifacts from naturalearthdata.com
earth    = NaturalEarth.naturalearth1("water")
borders  = NaturalEarth.borders()
airports = NaturalEarth.airports()
ports    = NaturalEarth.ports()

# initialize viewer with a coarse "raster"
earth |> Upscale(10, 5) |> viewer

# add other elements to the visualization
viz!(borders.geometry, color = "cyan")
viz!(airports.geometry, color = "black", pointsize=4, pointmarker='✈')
viz!(ports.geometry, color = "blue", pointsize=4, pointmarker='⚓')

# display current figure
Mke.current_figure()

image
image1200×900 102 KB

3

If I got it right, you want to make a figure with a 3D perspective of a topography and overlay on it a line network (not such a trivial job). From your words I think you made the first part, the 3D topography view but are struggling to overlay the lines. One workflow for doing this is:

  • Prepare the line network to be plotted with plot3. For this your lines must have a z dimension (it’s a 3D plot), and that z should match exactly the heights in the topography grid. This strep is accomplished by using grdtrack, which will interpolate the topo grid at the lines location.
  • Now do a grdview(...) followed by a plot3!(output_of_grdtrack, ...). The zscale and view angle are inherited from those set (automatically or not) in the grdview call.

Now the catch is that if the relief is very accentuated (likely not the case from your description), the lines lying on the invisible slopes will still be invisible because the plot3 programs has no hidden lines removal mechanism. The alternative here (2nd workflow) is to create first a flat image (grdimage + plot!) with the topography plus the line network and then drape this image on the 3D view using the draping option of grdview

This is an likewise example, but it uses 3D meshes instead of grids like you want.

4

Thanks for your answers. What I want is, I believe, simpler. Here is what I have so far:

using GMT
region = (16.30682249561994, 16.337816015500305, 48.19755410811093, 48.21355397447237)
surroundings = grdcut("@earth_relief_01s", region=region)
topographic_background = grdimage(
    surroundings,
    C=:earth,
    proj=:Winkel,
    colorbar=true,
)

This shows a 2D topographic map of a small part of a city (the resolution isn’t great but one thing at a time). I would “simply” like to have the street network displayed over this topographic background. I can find street network sources (either shape files or web services) but I have no clue about how to blend both in a way that looks good and that does not require me to draw the shape file “by hand” with street names etc.

5

All easy up to this point. Plotting shape files is easy. Just read them with

D = gmtread("path/to/shape.shp");

and add to the previous plot with

plot!(D)

Things start to give more work if you want to specify different line thickness for different streets. Then you have to split the D above by categories and use different line thicknesses. But street names must be plotted by hand meaning that positions and text orientation need to be specified. The plotting is easy (a call to text!) but preparing the data … not so much.

At this point you want to have a look at OpenStreetMapX to see if it helps.

6

Maybe take a look at Tyler

7

I thought in mentioning that too (but pointing to mosaic ofc :slight_smile:) but that will loose the relief topography info. Maybe a image blend of the OSM image and that topography can be a reasonable solution.

8

Made a quick test with the draping on 3D topography option. The SRTM topography is really too bumpy at this scale, but not obvious where to get a better one.

using GMT

region = (16.30682249561994, 16.337816015500305, 48.19755410811093, 48.21355397447237);
G = gmtread("@earth_relief_01s", region=region);
I = mosaic(G, provider=:OSM);
grdview(G, drape=I, view=(217,35), zsize=3, show=true)

bla
bla2295×1129 459 KB

and now with shading

grdview(G, drape=I, view=(217,35), zsize=3, show=true, shade=true)

bla_shade
bla_shade2295×1129 476 KB

9

Thanks a lot @joa-quim , this is definitely pointing me toward interesting directions!