Visualizing 3D Spatial Data With Pydeck
In this piece, I would like to introduce you to the Python package called PyDeck, which is a great tool to create 3D maps in Python. As an example, I will use a 3D building model data base covering the city of Budapest.
All images were created by the author.
Accessing the building data
The data is sourced from the publicly available Budapest Open Data Atlas, which you may access here. While the site is in Hungarian, once you find the section ‘Épületmagasság' and click on the button starting with ‘geojson', you should be able to download the geojson file, which is perfectly readable in English and with Python as well.
After a quick download, let's open the file and take a quick look:
# import geopandas
import geopandas as gpd
# parse the building height data set
gdf = gpd.read_file('epuletmagassag.geojson')
gdf.plot()
len(gdf)The output of this cell:
The image shows that the building footprints – more than 350k altogether – are scattered all around the city of Budapest. Now, let's take a quick look at the data structure more directly:
gdf.head(3)
For further use, two things to take from this table: the column DN containing the building height information, and the geometry containing the building footprint geometries.
Data preparation
Once you parse the GeoDataframe with 350k rows, you will also realize that for quick demonstrations, this is a bit too much. So, to make sure we can focus on a cool and efficient visual here, I decided to narrow down the data set to the central district of Budapest, District 5. For that, let's download the admin boundaries of this district using OSMnx and then carve out the corresponding footprints.
#Import all the libraries we use
import osmnx as ox
# Download the administrative boundary of the 1st district of Budapest
admin_district = ox.geocode_to_gdf('5th district, Budapest')
admin_district.plot()The outline of District 5:
# filter out the bulding footprints that fall into my selected districts
gdf2 = gpd.overlay(gdf, admin_district)
len(gdf2)This overlaid GeoDataFrame shows that we were successful in filtering the data – now we only have about 8k footprints.
By the way, the following quick histogram tells us that there are just three footprints above 60m in the district – which, after filtering the data frame, turn out to be two buildings – the Basilica and the Parliament. Matching this to my local knowledge, indeed, they are the two signature – and tallest – buildings in the area.
gdf2.DN.hist()
gdf2[gdf2.DN>60][['DN', 'name_1']]
Creating the visual
First, let's create a more basic version of the 3D building map by setting the coloring of each building to be red. With Pydeck, it's best to use RGB codes. Additionally, by changing the view_state, we can set up the center of the map as well as the zooming level.
import pydeck as pdk
layer = pdk.Layer(
'PolygonLayer',
data=gdf2,
get_polygon='geometry.coordinates',
extruded=True,
get_elevation='DN',
get_fill_color = [255, 0, 0],
)
view_state = pdk.ViewState(
latitude=47.499919,
longitude=19.0530236,
zoom=13,
)
# Create the 3D map
r = pdk.Deck(layers=[layer], initial_view_state=view_state)
r
Then, we can further fine-tune the plot, for instance, by changing the elevation scale to make the height differences more pronounced, deriving a scaling parameter, which we can incorporate in the coloring. This way, we can color the building models based on any parameter, or changing the line color as well.
import math
# do the scaling
max_val = max(gdf2.DN)
gdf2['scaled'] = [0.2 + (d/max_val) for d in gdf2.DN.to_list()]
# add the data
layer = pdk.Layer(
'PolygonLayer',
data=gdf2,
get_polygon='geometry.coordinates',
extruded=True,
get_elevation='DN',
elevation_scale=4,
get_line_color=[255, 255, 255],
get_fill_color="[255, 0, scaled * 255]"
)
# set the view state
view_state = pdk.ViewState(
latitude=47.499919,
longitude=19.0530236,
zoom=13,
)
# Create the 3D map
r = pdk.Deck(layers=[layer], initial_view_state=view_state)
r
Finally, tweak the view state a little further with camera direction settings, and save the output as an html file.
import math
# do the scaling
max_val = max(gdf2.DN)
gdf2['scaled'] = [0.2 + (d/max_val) for d in gdf2.DN.to_list()]
# add the data
layer = pdk.Layer(
'PolygonLayer',
data=gdf2,
get_polygon='geometry.coordinates',
extruded=True,
get_elevation='DN',
elevation_scale=4,
get_line_color=[255, 255, 255],
get_fill_color="[255, 0, scaled * 255]"
)
# set the view state
view_state = pdk.ViewState(
latitude=47.499919,
longitude=19.0530236,
zoom=14,
pitch=45, # controls the angle of the camera relative to the ground
bearing = 30 # controls the rotation of the camera around the vertical axis
)
# Create the 3D map
r = pdk.Deck(layers=[layer], initial_view_state=view_state)
r.to_html('test.html')
r
Conclusion
As this example shows, Python can be surprisingly useful in creating visuals that reveal various insights, such as building heights and other parameters, embedded in three dimensions.
If you would like to further advance your Geospatial Python skills, check out my brand new book, Geospatial Data Science Essentials – 101 Practical Python Tips and Tricks!
