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.
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')