this might be my finest work yet (fig. 1).

# method

overpass turbo^[1] for getting all the bakeries.

geopandas for clipping^[2] and reprojecting^[3].

seaborn jointplot for with kind='kde'

# all the code

import requests
import geopandas as gpd
import seaborn as sns
from shapely.geometry import box

overpass_url = "http://overpass-api.de/api/interpreter"

overpass_query = """
[out:json];
area["name"="France"]->.search_area;
(
  node(area.search_area)[shop=bakery];
);
out;
"""

# run the query
response = requests.get(overpass_url, 
                        params={'data': overpass_query})
data = response.json()

# grab lon lat from each element in list of dictionaries
lats = [d['lat'] for d in data['elements'] if d['type']=='node']
lons = [d['lon'] for d in data['elements'] if d['type']=='node']

# make geodataframe and clip to roughly mainland france
gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(lons, lats),
                       crs=4326)
gdf = gdf.loc[gdf.intersects(box(-5.5, 40, 8.5, 52))]

# reproject and add x, y coords
gdf = gdf.to_crs(gdf.estimate_utm_crs())
gdf['x'] = gdf.geometry.x
gdf['y'] = gdf.geometry.y

# do the plotting
jg = sns.jointplot(data=gdf,
                   x='x',
                   y='y',
                   kind='kde',
                   fill=False,
                   levels=20,
                   linewidths=0.5)

jg.ax_joint.set_aspect('equal')

# remove axes spines and labels
for ax in [jg.ax_joint, jg.ax_marg_x, jg.ax_marg_y]:
    ax.set_axis_off()

and yes. this is because croissant sounds a little bit like poisson, which is a distribution. i’m sorry.

# footnotes