Creating vector tiles for a road network
This short tutorial is meant as a follow up on this tutorial that explained the concept of vector tiles and how to create them in Python. In this previous article we already created a vector tiles, but we only looked at a simple case.
In this tutorial I will show you how to create vector tiles for a more advanced data set
As we saw there were 2 instances in which we had to go to the trouble of creating vector tiles.
1 The border of polygons are very detailed when compared to the size of the polygons. Take for example of contour lines. These can be very detailed and extremely long. Rendering an entire contour line in full detail is not efficient when we are just looking at a small section of it.
2 There are too many vectors and simply down sampling them is not good enough. This is the case for a road network for example. You do not which to just show a random sample of roads on low zoom level, instead you want to display the high ways or main roads.
This last case is more challenging then the first one. The first one can always be achieved by simplifying and tiling the vectors (in the way we did here). Case 2 however takes some creativity. How are we going to create low level views? Are we going to merge vectors? Display only certain types? Or a combination of the two?
In this tutorial we will have a look at case 2
The data set
We use an interesting data set containing the entire road network of the Netherlands (About 2 million on them!).
If you are interested or want to reproduce this example you can download the data set here. It is a large file of 1.5gb. If we read it as a geopandas
import geopandas as gpd
file = 'road_network'
sh = gpd.read_file(file)it even takes up 12Gb of RAM!
Now the question is… How do we visualize this data set in a meaningful way? Having our app download and display all segments would clearly be crazy. So we will need to apply our favorite trick vector tiles.
When exploring the data for a bit, you will likely discover that a road with 1 name can consist out of several segments. It would make sense to merge segments belonging to the same road. We could go further and merge roads with different names as well. The more we merge, the better we can simplify the view and the more performant we can render.
However, let’s say that I wish to visualize the roads based on their street name. This means that I cannot merge them further, since each vector with a different street name would need to be it’s own ‘unit’.
We can use the mergeGeometriesByColumn function from the ellipsis package. This function allows us to merge geometries in a geopandas that have the same value within a specified column.
sh = el.util.mergeGeometriesByColumn(features = sh, columnName = 'STT_NAAM', cores = 10)Here I use cores = 10 to parallelize the compute over 10 cores. (This data set is rather large so some optimization is needed to keep things feasible).
We now use street name, but depending on your use case other columns or derived columns will make more sense.
To construct meaningful lower zoom level views I split the data into highways, main roads and the rest. I will use this to show the highways first, and only when zooming in further the main roads and lastly the rest.
has_N = [ x != None and 'N' in x for x in sh['WEGNR_HMP'].values ]
sh_N = sh[has_N]
has_A = [ x != None and 'A' in x and not 'N' in x for x in sh['WEGNR_HMP'].values ]
sh_A = sh[has_A]
is_small = [ x == None for x in sh['WEGNR_HMP'].values ]
sh_small = sh[is_small]This gives us 3 geopandas dataframes. The high ways sh_A, the main roads sh_N and the rest sh_small.
Creating a layer
Before we continue manipulating the road network data we need to create a layer in Ellipsis Drive to have somewhere to write the vectors to. I can create a layer with timestamp within that layer using the following functions.
import ellipsis as el
token = el.account.logIn(MY_USERNAME, MY_PASSWORD)
pathId = el.path.add(pathType = 'vector', name = 'Styled roads server', token = token)['id']
timestampId = el.path.vector.timestamp.add(pathId = pathId, token=token)['id']Wit our layer created we can start adding our vectors. We start by adding the highways in the A_roads_merged geopandas. The highways need to be always visible. That means that we will add them to all zoom levels from 0 to the highest zoom we will use.
In this example I will use 13 as the highest zoom. If your vectors are more fine grained you can use a higher zoom.
As we add the vectors to lower and lower zoom levels I will use the util.simplifyGeometries function to simplify the geometries. On lower zoom we are further away so we can remove the unnecessary detail.
To construct a zoom level I always use the util.cutIntoTiles function. This functions splits up a vector if it intersects 2 or more tiles. This is important as vector tiles work best if all vectors are restricted to 1 tile alone.
With our simplify and cut to tiles function at the ready we can now execute the following code for the highways.
for zoom in np.arange(12, 14):
A_roads_merged_cut = el.util.cutIntoTiles(features=A_roads_merged, zoom=zoom, cores = 10, buffer = 0)
el.path.vector.timestamp.feature.add(pathId = pathId, timestampId = timestampId, token = token, features = A_roads_merged_cut, zoomLevels = [zoom], cores=10)
A_roads_merged = el.util.simplifyGeometries(A_roads_merged, tolerance= 100, cores = 10)
for zoom in np.arange(9, 12):
A_roads_merged_cut = el.util.cutIntoTiles(features=A_roads_merged, zoom=zoom, cores = 10, buffer = 0)
el.path.vector.timestamp.feature.add(pathId = pathId, timestampId = timestampId, token = token, features = A_roads_merged_cut, zoomLevels = [zoom], cores=10)
A_roads_merged = el.util.simplifyGeometries(A_roads_merged,tolerance = 10000, cores = 10)
for zoom in np.arange(0, 9):
A_roads_merged_cut = el.util.cutIntoTiles(features=A_roads_merged, zoom=zoom, cores = 10, buffer = 0)
el.path.vector.timestamp.feature.add(pathId = pathId, timestampId = timestampId, token = token, features = A_roads_merged_cut, zoomLevels = [zoom], cores =10)With this code done we added all highways to our layer. If we inspect it in the viewer we can see it was successful.
If we zoom in we can also observe that the detail of the roads increases.
Having the main roads done we will now add the main roads in N_roads_merged. Notice that we now add the vectors from zoom 13 till zoom 9, this is because we only want the main roads to show up if the user is zoomed in a bit further.
for zoom in np.arange(12, 14):
N_roads_merged_cut = el.util.cutIntoTiles(features=N_roads_merged, zoom=zoom, cores = 10, buffer = 0)
el.path.vector.timestamp.feature.add(pathId = pathId, timestampId = timestampId, token = token, features = N_roads_merged_cut, zoomLevels = [zoom], cores=10)
N_roads_merged = el.util.simplifyGeometries(N_roads_merged, tolerance=100, cores = 10)
for zoom in np.arange(9, 12):
N_roads_merged_cut = el.util.cutIntoTiles(features=N_roads_merged, zoom=zoom, cores = 10, buffer = 0)
el.path.vector.timestamp.feature.add(pathId = pathId, timestampId = timestampId, token = token, features = N_roads_merged_cut, zoomLevels = [zoom], cores=10)We can now see the secondary roads showing up in our map as well.
Now lastly we add the rest of the roads. There are a lot, and I mean a lot, of those so it is computationally very expensive. I use the cores parameter in the function to parallelize the operations, but then still it takes about 2 hours to complete this code.
We only add these smaller secondary roads to zoom 12 and 13.
zoom = 13
small_roads_merged_cut = el.util.cutIntoTiles(features=small_roads_merged, zoom=zoom, cores = 10, buffer = 0)
el.path.vector.timestamp.feature.add(pathId = pathId, timestampId = timestampId, token = token, features = small_roads_merged_cut, zoomLevels = [zoom], cores=10)
zoom = 12
small_roads_merged = el.util.simplifyGeometries(small_roads_merged, tolerance =100, cores = 10)
small_roads_merged_cut = el.util.cutIntoTiles(features=small_roads_merged, zoom=zoom, cores = 1, buffer = 0)
el.path.vector.timestamp.feature.add(pathId = pathId, timestampId = timestampId, token = token, features = small_roads_merged_cut, zoomLevels = [zoom], cores = 10)Final result
With all vectors added we now have a complete road map of the Netherlands. The further you zoom in the more detail you can see.
Conclusion
In this tutorial we created vector tiles using 2 techniques. We did not only simplify vectors but we also made a hierarchy ranking vectors in importance. This ranking allows us to make decisions on which vectors to add on lower zoom level and which vectors to skip.
The functions util.mergeGeometriesByColumn, util.simplifyGeometries and util.cutIntoTiles are 3 important helper functions to create vector tiles. They allow us to cut vectors on the tile borders, simplify them with a certain tolerance and to merge them base on properties of our choosing. These functions all have a cores parameter so that we can parallelize when needed.
