Calculating Building Density in R with OSM data
Introduction

OpenStreetMaps is a great source of spatial data. Most common programming languages have packages for downloading data from OSM.

In this tutorial I will show how to download building data using R’s osmdata package, perform density analysis and plot it using ggplot, and interactively using tmap.

This requires some knowledge of spatial data structures.

Getting the Data

The first step requires downloading data from the OSM api using osmdata. First, I will need to set a bounding box for OSM to search for data. Next, I will download administrative boundaries for making a choropleth map. Finally, I will download building data for the defined bounding box.

Setting a search area

I will be searching for data in the city of Izmir, Turkey. In order to define my search, I will use the getbb command. This command takes a location name as a string and returns a bounding box defined by OSM. getbb can also return a bounding box as a polygon, this is useful as I want to search for buildings as defined by administrative boundaries.

Because the function returns a matrix, you can use your own polygon to define the OSM opq lookup (below), but it is important to make sure the CRS of the data is EPSG:3857.

Getting the bounding polygon

Retrieve boundary lines

In this project I am looking for neighborhood boundaries of districts in Izmir. Because administrative level boundaries are not standardized across countries, it may take some research to figure out what administrative levels are needed. Luckily, there his a helpful tool that lets you find OSM boundaries within an interactive map.

The numbers next to territory names signify the administrative level that is saved within OSM. In my case, I am looking for boundaries at level 8. I can pass this information using the osmdata request feature, which is a wrapper for the OSM API.

Retrieve OSM boundary lines

(Note: OSM has many more map features available and can be queried using the right key:value pairs.)

The function returns an osmdata object, which is a list of sf dataframes. Each dataframe is composed of a different geometry type, this means that returned objects can be large. In order to get the right data, you need to know which geometry you are after.

In this case, the shape I need are multipolygon geometry types. OSM querys can end up returning large objects. In order to keep some space free, I will drop the extra dataframes. I’ll also clean up some of the data.

Cleaning boundary data

Now I want to check to see if I have the right boundaries.

OSM Boundary lines for original bounding polygon

Looks good, but it has also returned the boundaries for all the outlying districts, as they probably fall along the bounding polygon line. That is something I can deal with later.

Retrieve Buildings

Using the same process as above, query OSM for all buildings within the original bounding polygon: iz_bbox . The building polygons object will be about 100x larger than my boundaries object, so I’ll get rid of the excess data again. This time I will keep the polygon sf object.

Then I’ll plot the result.

Plotting building polygons over boundary polygons

The buildings returned don’t really follow the borders of any polygon, so now I’ll do some spatial trimming.

Clipping point data

In order to filter out all of the outlying buildings, I need to calculate them as points, which can be done by finding their centroids. I’ll create a new dataframe with the building centroids. Before I do that, I want to calculate the area of all the buildings, so the data carries over with the newly created dataframe. I won’t need the building polygons from this point on.

I will also need to create a polygon geometry in order to filter the outlying buildings. I can do this using the original bounding box iz_bbox . But, because this object is a matrix, I need to convert it to an sf object. This can be done using sf_headers::sf_polygon .

Polygon in red is the target area, I’ll can use this polygon for clipping building data

Now I have a polygon that I can filter points with. First, I will have to create a dataframe that will allow me to filter. setting st_join join variable to st_within will give me a column indicating whether a point falls within the polygon dataframe. The column id is created, with a 1 indicating that a point falls within the polygon. I’ll plot to check the result.

Building centroids as points within the bounding polygon

It worked, so now I’ll filter the right points, and then create a new dataframe from the iz_polys dataframe, but with the filtered points.

Using the more granular iz_polys dataframe, I can aggregate and sum the total area of all buildings within each district. I will create a final dataframe that contains the total building area as well as the total area of the district polygon.

Now I’ll create a nice interactive map using tmap .

Screenshot of interactive map. See the real thing here: https://greggsaldutti.netlify.app/post/2021-05-05-calculating-building-density-in-r-with-osm-data/

Conclusion

So here we have an areal density map of buildings as a percentage total land area. R is a great tool when you have lots of data to wrangle, or need to access data programatically. Other programs, such as QGIS, may be more user friendly, but R gives you much more control, and makes great looking maps too.

Notes: This post is based off of Matt Herman’s post about trees in NYC, and this great open-source GIS course using Python.
['Gregg Saldutti']
['Rstats', 'Openstreetmap', 'Spatial Analysis', 'Mapping', 'GIS']