In my last post I prepared a geopackage with all peaks of the Alps that can be found on Openstreetmap. But how can I add columns about countries or the corresponding mountain areas? Well, a peak can be on a border and belong to more than one country, but it should only be within one mountain area. This means we need two different approaches. Of course, it is easy to adapt the code to use with any other kind of polygon (region, city, etc.).
The following code is available as Jupyter notebook on GitHub.
The data used in this notebook is © OpenStreetMap contributors under the Open Database License.
Add countries
For this, I downloaded the countries of the area with QuickOSM (key: admin_level, value: 2).
Since a peak can belong to several countries, I want one column per country with True/False. Note that this approach is slow.
import pandas as pd
import geopandas as gpd
import os
folder = 'alpinepeaks/'
# Load the peaks geopackage
gdf = gpd.read_file(folder + "peaks.geojson")
gdf.set_index('full_id', inplace=True)
# Get rid of the double Monaco
countries.replace('Monaco (territorial waters)', 'Monaco', inplace=True)
countries = countries.dissolve(by='name:en')
# The last line also sets 'name:en' as index
# That means we only need index + geometry
countries = countries[[ 'geometry']]
# Add a buffer of 50 m to countries
# (first reproject from WGS84 to UTM)
# to make sure that peaks on the border really get both countries.
countries = countries.to_crs("EPSG:32634")
countries.geometry = countries.geometry.buffer(50)
# Back to WGS84
countries = countries.to_crs("EPSG:4326")
# I create one column per country with True or False. The loop is slow!
for c in countries.index:
print('processing', c)
gdf[c] = gdf['geometry'].within(countries.loc[c].geometry)
# Save
gdf.to_file(folder + "peaks2.geojson", driver='GeoJSON')
Add mountain area
I assume that a peak is only in one mountain area. This is much faster, as I only add one column.
From openstreetmap I downloaded relations with region:type=mountain_area.
# Open mountain regions geopacke
area = gpd.read_file(folder + 'mountain_region.gpkg')
# I define a function so I don't need to loop over the data frame.
def mountain_area(point):
for a in area.index:
if point.within(area.loc[a].geometry):
return area.loc[a]['name']
return None
# Map the function
gdf['mountain_area'] = gdf['geometry'].map(mountain_area)
# Save
gdf.to_file(folder + "peaks2.geojson", driver='GeoJSON')