UFO Reports Across the United States

Mapping the hotspots and modeling the rhythm of when people look up

This notebook is a portfolio case study by AbleVLabs. It works end to end with a real dataset: roughly 80,000 geolocated UFO reports collected by the National UFO Reporting Center (NUFORC), accessed through the public TidyTuesday mirror.

One honest note up front. This project models the pattern of human reports. It does not, and cannot, predict real anomalous events. Every trend and every spike in this data has an ordinary explanation: when the reporting website existed, when people were outside on warm nights, and which holidays put crowds under the sky. That framing is the whole point. Good data science names what the numbers actually measure.

1. The question

Three questions drive the analysis:

1. Where do reports cluster across the country?
2. When do they happen, across the season, the week, the clock, and the calendar?
3. How has the volume of reports changed over time, and what does a forecast of that volume actually represent?

The deliverable is a set of clear patterns plus a small, transparent model that turns a month, a weekday, and an hour into a relative likelihood of a report.

2. The data

Source. NUFORC (nuforc.org), the longest running public UFO report collector in the United States, mirrored by the TidyTuesday project (rfordatascience/tidytuesday, 2019-06-25 release).

Shape. The extract carries about 80,000 US reports spanning 1925 to 2014, with a timestamp, a city and state, latitude and longitude, a reported shape, an encounter length, and a free-text description.

The full extract was aggregated at the source into the distributions used below, and a random 900-report geolocated sample was drawn for the map. Both the aggregates and the sample are the real NUFORC records, not synthetic stand-ins.

import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

# dark theme to match the AbleVLabs site
plt.rcParams.update({
    'figure.facecolor':'#0d0d12','axes.facecolor':'#12121a',
    'savefig.facecolor':'#0d0d12','text.color':'#e8e8ef',
    'axes.labelcolor':'#e8e8ef','xtick.color':'#a9a9c0','ytick.color':'#a9a9c0',
    'axes.edgecolor':'#2a2a3a','grid.color':'#20202c','font.size':11,
    'axes.titlecolor':'#e8e8ef','figure.dpi':110,
})
CY, PU, GR, AM = '#22d3ee', '#a855f7', '#34d399', '#f59e0b'

D = json.load(open('ufo.json'))
meta = D['meta']
print('source :', meta['source'])
print('span   :', meta['year_min'], 'to', meta['year_max'])
print('raw US :', f"{meta['us_total']:,} reports")

source : NUFORC (National UFO Reporting Center) via TidyTuesday mirror
span   : 1925 to 2014
raw US : 65,111 reports

3. First look: the geolocated sample

The 900-report sample loads into a tidy DataFrame. Each row is one real report with a location, a year, an hour of the day, and a reported shape. This is the frame used for the row-level work; the full-population counts are used for the big-picture distributions.

df = pd.DataFrame(D['sample']).rename(columns={
    'la':'lat','lo':'lon','st':'state','y':'year','h':'hour','sh':'shape'})
print(df.shape)
df.head(8)

(900, 6)

	lat	lon	state	year	hour	shape
0	39.96	-83.00	OH	2001	20	light
1	43.97	-75.91	NY	2002	21	other
2	36.06	-94.16	AR	2013	21	chevron
3	41.80	-71.89	CT	2005	8	oval
4	36.40	-93.74	AR	1999	21	light
5	42.87	-106.31	WY	2011	23	light
6	44.36	-98.21	SD	2009	0	light
7	39.56	-76.07	MD	2009	23	circle

# a quick profile of the sample
display(df[['lat','lon','year','hour']].describe().round(2))
print('\nstates represented :', df['state'].nunique())
print('distinct shapes    :', df['shape'].nunique())

	lat	lon	year	hour
count	900.00	900.00	900.00	900.00
mean	38.56	-94.92	2003.97	15.52
std	5.45	17.67	10.95	7.81
min	21.42	-158.19	1950.00	0.00
25%	34.16	-112.07	2001.00	10.00
50%	39.22	-89.52	2007.00	19.00
75%	42.13	-80.32	2011.00	21.00
max	64.84	-68.84	2014.00	24.00

states represented : 50
distinct shapes    : 21

4. Cleaning

Real report logs are messy. Before any analysis the source pass kept only rows that sit in a real US state and carry a usable location, and it standardized the timestamp into a year, a month, an hour, and a weekday. The table below is the audit trail: what came in, what was set aside, and what remained to analyze.

cl = D['clean']
audit = pd.DataFrame({
    'stage':['raw rows returned','dropped: non-US or missing geo','US analyzable reports'],
    'rows':[cl['raw'], cl['non_us_or_missing_geo'], cl['us_analyzable']],
})
audit['pct_of_raw'] = (audit['rows']/cl['raw']*100).round(1)
audit

	stage	rows	pct_of_raw
0	raw rows returned	80327	100.0
1	dropped: non-US or missing geo	15216	18.9
2	US analyzable reports	65111	81.1

About 81 percent of the returned rows resolve cleanly to a US state with coordinates. The rest are international reports or rows with a missing or unparseable location, and they are set aside so the geography is honest.

5. Where: the hotspots

Two views of geography. First the raw leaderboard of states, then a scatter of the real 900-report sample on the map. The coasts and the Sun Belt dominate, which tracks population and clear-sky nights more than anything otherworldly.

states = pd.Series(D['byState']).sort_values(ascending=False)
top = states.head(15)
fig, ax = plt.subplots(figsize=(9,4.6))
ax.bar(top.index, top.values, color=CY, edgecolor='#0d0d12')
ax.set_title('Top 15 states by report count (full NUFORC US extract)')
ax.set_ylabel('reports')
ax.yaxis.set_major_formatter(mticker.StrMethodFormatter('{x:,.0f}'))
for i,v in enumerate(top.values):
    ax.text(i, v+80, f'{v:,}', ha='center', va='bottom', fontsize=8, color='#a9a9c0')
plt.tight_layout(); plt.show()

# scatter map of the real sample over the continental US
m = df[(df.lon>-125)&(df.lon<-66)&(df.lat>24)&(df.lat<50)]
fig, ax = plt.subplots(figsize=(9,5.4))
ax.scatter(m.lon, m.lat, s=16, c=CY, alpha=0.55, edgecolors='none')
ax.set_title(f'Reported UFO locations, random sample of {len(m)} real reports')
ax.set_xlabel('longitude'); ax.set_ylabel('latitude')
ax.set_xlim(-125,-66); ax.set_ylim(24,50)
ax.grid(True, alpha=0.25)
plt.tight_layout(); plt.show()

The outline of the country appears without any map layer underneath it. That is the signature of a population-driven process: reports trace where people live, the coasts, the Great Lakes, Texas, and Florida.

6. When: the season

Aggregating every report by calendar month exposes a strong summer bulge. July is the single busiest month. Warm nights, vacations, fireworks, and time spent outdoors all push the count up.

mo = pd.Series({int(k):v for k,v in D['byMonth'].items()}).sort_index()
names = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
fig, ax = plt.subplots(figsize=(9,4))
bars = ax.bar(names, mo.values, color=PU, edgecolor='#0d0d12')
bars[mo.values.argmax()].set_color(AM)
ax.set_title('Reports by month (summer peaks, winter dips)')
ax.set_ylabel('reports')
ax.yaxis.set_major_formatter(mticker.StrMethodFormatter('{x:,.0f}'))
plt.tight_layout(); plt.show()
print('July vs February:', round(mo[7]/mo[2],2), 'x')

July vs February: 2.04 x

7. When: the clock

The hour-of-day distribution is the sharpest pattern in the whole dataset. Reports climb through the evening and peak at 9 to 10 pm, then fall through the small hours. This is the "after dinner, before bed, still dark, still awake" window.

hr = pd.Series({int(k):v for k,v in D['byHour'].items()}).sort_index()
fig, ax = plt.subplots(figsize=(9,4))
cols = [AM if h in (21,22) else GR for h in hr.index]
ax.bar(hr.index, hr.values, color=cols, edgecolor='#0d0d12')
ax.set_title('Reports by hour of day (9 to 10 pm is the peak)')
ax.set_xlabel('hour (24h)'); ax.set_ylabel('reports')
ax.set_xticks(range(0,24,2))
ax.yaxis.set_major_formatter(mticker.StrMethodFormatter('{x:,.0f}'))
plt.tight_layout(); plt.show()
print('9pm vs 8am:', round(hr[21]/hr[8],1), 'x')

9pm vs 8am: 14.3 x

8. When: the week and the calendar

Saturday leads the week, with Friday and Sunday close behind: the leisure pattern again. And a single calendar day towers over all others. The Fourth of July is the busiest date in the dataset by a wide margin, followed by other round or holiday dates that also happen to draw crowds outdoors (and fill the sky with fireworks and lanterns).

wd = pd.Series({int(k):v for k,v in D['byWeekday'].items()}).sort_index()
wdn = D['model']['wd_names']
fig, axes = plt.subplots(1,2, figsize=(11,4))
b = axes[0].bar(wdn, wd.values, color=CY, edgecolor='#0d0d12')
b[wd.values.argmax()].set_color(AM)
axes[0].set_title('Reports by weekday'); axes[0].set_ylabel('reports')
axes[0].yaxis.set_major_formatter(mticker.StrMethodFormatter('{x:,.0f}'))

days = D['topDays'][:8]
lbl = [d[0] for d in days][::-1]; val=[d[1] for d in days][::-1]
cc = [AM if l=='7-4' else PU for l in lbl]
axes[1].barh(lbl, val, color=cc, edgecolor='#0d0d12')
axes[1].set_title('Busiest calendar dates (month-day)'); axes[1].set_xlabel('reports')
plt.tight_layout(); plt.show()
print('Fourth of July reports:', dict(D['topDays'])['7-4'])

Fourth of July reports: 1208

9. The trend: modeling report volume over time

Annual report volume rises steeply, especially from the mid-1990s on. The critical insight is what this curve measures. NUFORC's web form arrived in the internet era, so the climb is largely a story of access and awareness, not of more objects in the sky.

A log-linear model is fit to the modern reporting regime (1995 to 2013). The 2014 slice is a partial year and is held out of the fit. The model is then backtested by training on 1995 to 2010 and predicting the held-out 2011 to 2013.

tr = D['trend']
yrs = np.array(tr['fit_years']); act = np.array(tr['fit_actual'])
# log-linear fit reproduced here for transparency
b1,b0 = np.polyfit(yrs, np.log(act), 1)
fit = np.exp(b0 + b1*yrs)
growth = (np.exp(b1)-1)*100
print(f'implied growth: {growth:.1f}% reports per year')

bt = tr['backtest']
print(f"backtest {bt['test_years'][0]}-{bt['test_years'][-1]}:",
      f"MAE={bt['mae']:,} reports, MAPE={bt['mape']}%")
print('actual :', bt['actual'])
print('predict:', bt['pred'])

implied growth: 10.2% reports per year
backtest 2011-2013: MAE=1,081 reports, MAPE=18.4%
actual : [5107, 7357, 7037]
predict: [6810, 7554, 8379]

fig, ax = plt.subplots(figsize=(9.4,4.8))
ally = np.array(tr['all_years']); alla = np.array(tr['all_actual'])
ax.plot(ally, alla, 'o-', color=CY, ms=4, lw=1.4, label='actual reports')
ax.plot(yrs, fit, '--', color=AM, lw=2, label='log-linear fit (1995-2013)')
fy = tr['forecast_years']; fp = tr['forecast_pred']
ax.plot(fy, fp, 's--', color=PU, ms=5, lw=1.6, label='forecast 2014-2018')
ax.set_title('Annual report volume: actual, fit, and forecast')
ax.set_xlabel('year'); ax.set_ylabel('reports')
ax.yaxis.set_major_formatter(mticker.StrMethodFormatter('{x:,.0f}'))
ax.legend(facecolor='#12121a', edgecolor='#2a2a3a', labelcolor='#e8e8ef')
ax.grid(True, alpha=0.25)
plt.tight_layout(); plt.show()

The forecast projects report volume forward at roughly 10 percent per year. It is a projection of reporting behavior, contingent on NUFORC staying the go-to venue and on public interest holding steady. It says nothing about physical events. A backtest error near 18 percent is reasonable for a single-trend model on noisy annual counts.

10. A likelihood model for "when"

The seasonal, weekly, and hourly patterns combine into a simple, transparent model. Each dimension contributes a factor equal to its observed share divided by a uniform share, so a factor of 1.0 means "an average slice" and 2.0 means "twice as likely as random." The month factor divides out the differing month lengths so it is a fair per-day rate.

Multiplying the three factors gives a relative likelihood index for any month, weekday, and hour. This is a naive-independence model, stated plainly, and it is exactly what the interactive predictor on the dashboard uses.

M = D['model']
mf = {int(k):v for k,v in M['monthFactor'].items()}
wf = {int(k):v for k,v in M['weekdayFactor'].items()}
hf = {int(k):v for k,v in M['hourFactor'].items()}

def likelihood_index(month, weekday, hour):
    """Relative likelihood of a report vs an average random hour (1.0 = average)."""
    return round(mf[month]*wf[weekday]*hf[hour], 2)

# the single most likely slot in the whole model
pk = M['peak']
print(f"peak slot -> month {pk['month']}, {M['wd_names'][pk['weekday']]}, {pk['hour']}:00")
print(f"peak index: {pk['index']}x baseline")
print()
print('July, Saturday, 9pm  :', likelihood_index(7,6,21), 'x')
print('Feb,  Tuesday,  8am  :', likelihood_index(2,2,8), 'x')
print('Oct,  Friday,  10pm  :', likelihood_index(10,5,22),'x')

peak slot -> month 7, Sat, 21:00
peak index: 5.9x baseline

July, Saturday, 9pm  : 5.9 x
Feb,  Tuesday,  8am  : 0.17 x
Oct,  Friday,  10pm  : 3.59 x

# hour x month likelihood heatmap (weekday held at Saturday)
sat = 6
grid = np.array([[mf[mo]*wf[sat]*hf[h] for mo in range(1,13)] for h in range(24)])
fig, ax = plt.subplots(figsize=(9.2,5))
im = ax.imshow(grid, aspect='auto', origin='lower', cmap='inferno')
ax.set_xticks(range(12)); ax.set_xticklabels(names)
ax.set_yticks(range(0,24,2)); ax.set_yticklabels(range(0,24,2))
ax.set_xlabel('month'); ax.set_ylabel('hour (24h)')
ax.set_title('Relative likelihood index on a Saturday (brighter = more likely)')
cb = fig.colorbar(im, ax=ax); cb.set_label('index (1.0 = average)')
plt.tight_layout(); plt.show()

The heatmap says the quiet: a February morning is almost dead. And the loud: a July evening around 9 to 10 pm is where reports concentrate. The single most likely slot the model can name is a July Saturday at 9 pm, and the busiest date on the real calendar is the Fourth of July. If the question is "when is the next likely wave of reports," the honest answer the data supports is: the next warm-weekend holiday evening.

11. Conclusions and limitations

What holds up.

• Geography follows population. California, Washington, Florida, Texas, and New York lead.
• The clock is the strongest signal: reports peak at 9 to 10 pm.
• The season peaks in July, the week peaks on Saturday, and the Fourth of July is the single busiest date by a wide margin.
• Report volume grew about 10 percent per year through the online era.

What this is not.

• It is not evidence about real objects. Every pattern matches ordinary human behavior: internet access, leisure time, warm nights, and holidays.
• The trend forecast projects reporting behavior, not events, and assumes the reporting channel stays stable.
• The likelihood model assumes the three time dimensions act independently, which is a simplification.

Why it still matters. The same toolkit, geospatial aggregation, temporal decomposition, a backtested trend model, and a transparent likelihood score, is exactly what powers demand forecasting, staffing models, and anomaly detection in industry. The subject is playful. The method is real.

Data: NUFORC via the TidyTuesday mirror. Analysis and visuals by AbleVLabs.