From Simple Points to Angry Red Blobs: Kernel Density Mapping in Python
Let's say that there are 3,023 squirrels in Central Park. There are probably way more, but those are the ones that someone saw and recorded the location. If we put them on a map, we get 3,023 dots. This is accurate, but still a bit.
A scatter plot can tell us where individual observations happened, but it struggles with the question people usually care about:
Where are the squirrels actually concentrated?
That is where kernel density estimation, or KDE, enters the picture.
In this post, we will take one real dataset from raw point coordinates to a continuous density surface and finally to a classified hotspot map. Along the way, we will see why coordinate systems matter, why bandwidth controls the entire story, and why two Python libraries can produce slightly different answers from the same squirrels.
The full journey is:
points → kernels → overlapping influence → density surface → hotspots
We will not need any advanced GIS packaes. Just Python together with a few plotting libraries, and a lot of squirrel data.
Meet the Dataset: The 2018 Central Park Squirrel Census
In October 2018, more than 300 volunteers spread across Central Park and recorded Eastern gray squirrel sightings. The resulting dataset contains 3,023 observations, including coordinates, fur colour, activities, and several behavioural notes.
It is a wonderfully specific dataset, with no sensible real world usage. But, it is also well suited to a KDE demonstration.
We can load the public copy directly with Pandas:
For completeness we make a local copy and use that if it already exists. The dataset contains plenty of entertaining fields. If you are interested, there are
long
lat
unique_squirrel_id
hectare
shift
date
hectare_squirrel_number
age
primary_fur_color
highlight_fur_color
combination_of_primary_and_highlight_color
color_notes
location
above_ground_sighter_measurement
specific_location
running
chasing
climbing
eating
foraging
other_activities
kuks
quaas
moans
tail_flags
tail_twitches
approaches
indifferent
runs_from
other_interactions
lat_long
zip_codes
community_districts
borough_boundaries
city_council_districts
police_precincts
But we only need two of them for now:
lat
long
Those two columns describe individual locations. Our goal is to turn those isolated locations into an estimate of spatial concentration.
Before we can do that, however, we need to fix our coordinate system.
Latitude and Longitude to Metres
A KDE bandwidth needs a meaningful distance. If we say that each sighting influences an area within roughly 100 metres, then the number 100 should actually mean 100 metres. Raw longitude and latitude values do not give us that. One degree of longitude also covers a different physical distance depending on latitude.
For a compact study area such as Central Park, a local metre-based approximation is sufficient:
This centres the coordinates on the dataset and converts them into local x and y distances measured in metres.
It is not a replacement for a proper projected coordinate reference system in a large production GIS workflow. But for a small area and an explanatory example, however, it gives us exactly what we need: distances that behave like distances.
That means a bandwidth of 60 now means roughly 60 metres, rather than “0.0006 degrees and good luck explaining that to anyone.”
Giving the Squirrels Somewhere to Stand
A density surface without geographic context is just a colourful rectangle.
To keep the maps readable, we will use a simplified outline of Central Park and represent the Jacqueline Kennedy Onassis Reservoir as a rotated ellipse. This is intentionally lightweight. The goal is to explain KDE, not disappear into a side quest involving shapefiles, API keys, invalid polygons, and one multipolygon that refuses to cooperate.
The following helpers create the common map background used throughout the post:
The park boundary is only an approximation, but it provides enough visual context to distinguish “a hotspot in the southern lawns” from “a hotspot somewhere inside a rectangle.” Now we can look at the raw data.
The Honest but Unhelpful Scatter Plot
Plotting all 3,023 sightings gives us the most literal view of the census:
Raw point data over the real park outline gives us then the following scatter plot.
The plot is truthful. But really, it is also a brown cloud shaped vaguely like Central Park.
We can see that some areas contain more sightings than others, but the map does not tell us where one cluster ends, where another begins, or how strongly each part of the park is supported by nearby observations.
This is the point-cloud problem, namely the information is present, but our eyes have to perform the analysis manually.KDE replaces that visual guesswork with a repeatable calculation.
The Kernel: Every Squirrel Gets a Tiny Fog Machine
The basic idea behind KDE is simple. Instead of treating a sighting as an infinitely small point, we place a smooth bump over it. The bump is highest at the observation and fades with distance. A common choice for this would be a Gaussian curve.
The width of that bump is controlled by the bandwidth, usually written as h.
Here is the same Gaussian kernel with two different bandwidths:
As a plot this gives us the following two little bumps
A small bandwidth creates a narrow, sharp bump. In this case, the sighting only influences its immediate surroundings. A large bandwidth on the other hand creates a wide, gentle bump. In this case, the same sighting now influences a much larger part of the park.
The observation has not changed, rather its influence that we associate with it has.
That is the most important idea in this case:
Bandwidth is not a cosmetic setting. It determines the spatial story the map tells.
A 60-metre bandwidth might distinguish separate benches. A 220-metre bandwidth might merge several lawns into one broad zone. Neither value is automatically correct. The right choice depends on the scale of the process we are trying to understand.
Where the Fog Overlaps, Density Builds
Of course, one kernel is not very exciting, but at three thousand kernels this might change a bit.
Whenever several sightings occur close together, their Gaussian bumps overlap. The individual contributions are added, producing a higher total density in that area.
We can illustrate the idea with visible circles around a few sample observations:
In this plot the little red circle around 70 observations indicate their influence, on the left plot the influence is chosen as 60m, while in the right plot it is chosen as 120m.
The circles are only a diagram, not the actual Gaussian calculation, but they reveal the core mechanism. A patch of ground covered by several circles receives several contributions. A patch covered by none receives almost nothing. KDE performs this summation smoothly across a fine grid. It is essentially asking the same question at thousands of locations:
How much support does this cell receive from all nearby sightings?
Once every cell has an answer, we have moved from a point cloud to a coloured canvas.
Building the Density Surface with SciPy
SciPy provides gaussian_kde, which can estimate a smooth density from the point coordinates. First, we create a regular grid covering the park. Then we evaluate the KDE at every grid location:
There is one technical wrinkle hiding inside this helper. The method scipy.stats.gaussian_kde does not treat the bandwidth as a simple circular radius. It starts with the covariance structure of the data and applies a scalar factor to that covariance matrix. This essentially means that it looks at how the points vary in both dimensions and scales everything accordingly.
In our example, Central Park is long and narrow. Hence, the sighting cloud follows that shape. As a result, SciPy’s kernels can stretch along the park’s main axis rather than remaining perfectly circular.
The conversion in scipy_kde_surface() gives us a repeatable metre-based control, but it does not turn SciPy’s estimator into an isotropic spatial kernel, i.e. one where the influence would be a perfect circle. That distinction will matter later.
For now, let us see what happens when we change the bandwidth.
One Dataset, Two Completely Different Maps
We will calculate two surfaces:
a tight 60-metre version;
a broad 120-metre version.
Both use the same sightings, grid, and plotting code:
This gives us two very different types of surfaces/
Note that the observations are exactly the same, the only difference is our choice of influence.
The 60-metre surface preserves many small local peaks. It is detailed, but some of that detail may reflect survey paths, volunteers pausing in one place, or random clusters of nearby observations.
The 120-metre surface suppresses those small fluctuations. It reveals broad zones, but it also merges locations that may be meaningfully different on the ground.
This is the classic KDE trade-off:
too little smoothing gives us noise wearing a colourful hat;
too much smoothing gives us one enormous blob with excellent self-confidence.
For the rest of the post, we will use a middle value of 90 metres. It retains the park’s larger structure without assigning a separate hotspot to every bench, path junction, and suspiciously squirrel-shaped rubbish bin.
The Continuous Surface
With the 90-metre bandwidth selected, we can produce the main density surface:
This is a middle ground between our two previous choices.
This is now a continuous raster-like surface. Every grid cell has a density estimate, including cells where no squirrel was directly observed. That does not mean KDE has invented new sightings. It means each cell receives support from observations nearby.
A high value says:
Many sightings occurred close to this location.
A low value says:
Few sightings occurred nearby, or the nearest sightings were relatively far away.
The result is much easier to interpret than the raw scatter plot, but it is still a continuous gradient. If we want clear hotspot classes, we need one more step.
Turning the Surface into Hotspots
A smooth density map is useful for exploration. A classified map is often easier to communicate.
We can divide the surface using percentile thresholds. In this example:
values above the 75th percentile become a low hotspot tier;
values above the 90th percentile become a stronger tier;
values above the 95th percentile become the highest tier.
The code below calculates those thresholds and assigns each grid cell to a class:
This gives us only 4 different zones.
Now the map expresses a much clearer message. The darkest background represents lower-density areas. The yellow, orange, and red regions show increasingly concentrated clusters of sightings.
This directly tells you where the main sightings happened. But there is an important limit:
A percentile is a ranking, not a significance test.
The 95th percentile means that a cell belongs to the highest five per cent of values on this particular surface. It does not prove that the cluster is statistically unusual, biologically important, or impossible under a random spatial process.
For formal hotspot testing, we would need a spatial statistic such as Getis-Ord Gi*. PySAL’s esda package supports that kind of analysis.
For a practical first look at where observations concentrate, percentile classes are simple, transparent, and useful. They just should not be promoted to scientific certainty because they look convincing.
One Kernel, Two Libraries
SciPy is not the only Python option for KDE. Scikit-learn provides sklearn.neighbors.KernelDensity, which has a particularly convenient property for spatial work: its bandwidth is expressed directly in the units of the input coordinates.
Because our coordinates are measured in metres, bandwidth=90 means a circular 90-metre Gaussian kernel.
Let us compare the two implementations:
Looking at the two surfaces we see clear differences
The surfaces have some similarity in where there are hotspots and no sightings, but they are far from being identical.
SciPy’s gaussian_kde uses the covariance structure of the whole point cloud. Since Central Park is elongated, its smoothing stretches in that direction. Which means that observations along the length of Central Park have more interactions and our areas broadly follow the same stretched outline.
Scikit-learn’s KernelDensity on the other hand, uses the explicit isotropic bandwidth supplied to it. The smoothing radius is circular in the coordinate space, which preserves more local spatial detail. Which shows up as hot spots that are much more local and less shaped like the park itself.
The practical distinction is useful:
Library | Best suited to |
scipy.stats.gaussian_kde | Quick exploratory density estimation with covariance-aware smoothing |
sklearn.neighbors.KernelDensity | Explicit spatial bandwidths, alternative kernel shapes, and integration with scikit-learn workflows |
Neither library is universally better. The more important question is whether the estimator’s assumptions match the map you think you are making.
Where KDE Can Mislead You
KDE is simple to run, which makes it dangerously easy to trust. Before turning the result into a report, dashboard, or some red blob for a management presentation, there are several things worth checking.
Sampling Bias Does Not Disappear
KDE smooths the data you give it. But it does not know how the data was collected. The Squirrel Census used a coordinated survey design, which makes it unusually suitable for this example. Many real datasets are less balanced.
A hotspot in app-reported wildlife observations might simply be an area with more app users. A hotspot in reported incidents might reflect easier access to reporting services. A hotspot in delivery requests might mirror population density rather than unusual demand.
KDE cannot separate “more events happened here” from “we observed this area more carefully.”
Bandwidth Must Be Reported
A hotspot map without a stated bandwidth hides its most important modelling decision.
Whenever possible, express the bandwidth in meaningful physical units and show multiple alternatives. A side-by-side comparison often communicates the uncertainty more honestly than one polished final map.
Percentiles Are Not Statistical Tests
Percentile classes are descriptive. They tell us which cells rank highest within the estimated surface.
They do not tell us whether those cells are statistically significant under a null model. For completeness here is the hotspot version of the scikit-learn KDE.
Boundaries and Water Still Exist
A Gaussian kernel has no knowledge of fences, lakes, roads, cliffs, or administrative borders.
Near the edge of a study area, kernels extend beyond the boundary, which can reduce estimated density inside it. This is the classic edge-effect problem.
The reservoir creates another obvious example. A smooth surface may spread density across water because the algorithm sees coordinates, not squirrel swimming preferences.
More advanced workflows can apply masks, boundary corrections, network distances, or habitat constraints. The basic KDE shown here does none of those things.
Final Thoughts
We started with 3,023 individual points and ended with a map that describes spatial concentration. The calculation itself is not the difficult part. The code required to fit a KDE and evaluate it on a grid is compact.
The real decisions were:
how to represent distance;
what bandwidth to use;
whether the kernel should be isotropic;
how to classify the continuous surface;
how to account for sampling, edges, and geographic barriers.
Those choices determine whether a KDE map becomes a useful analytical tool or merely an angry red blob with a colourbar.
The same workflow applies far beyond squirrels. Crime reports, disease cases, animal sightings, service calls, traffic incidents, deliveries, and sensor events all begin as point locations. KDE provides a way to move from individual events to a broader spatial pattern.
We just have to remember that the smooth surface is an estimate, not a discovery machine.
It can show where the data are concentrated. It cannot tell you why the squirrels chose that lawn.
Comments