Геолокация острова по фото: геометрия, OpenStreetMap и CUDA

16-08-2026

Автор материала подчёркивает: работа выполнена вручную, без использования LLM.

Разбор посвящён решению OSINT-задачи gralhix 004, подготовленной Sofia Santos | Gralhix.

Все файлы кода и итоговый отчёт с инструкциями доступны на GitHub.

Условие задачи

main

Дано фото курорта, расположенного на острове.

a) Как называется курорт?
b) Каковы координаты острова?
c) В каком направлении была направлена камера в момент съёмки?

Решать задачу через Google Lens означало бы упустить интересную возможность, поэтому был выбран путь через математику и программирование.

a] Метаданные

Первым делом проверены метаданные файла — команда запущена в Linux Void:

> exiftool main.png

File Type                       : WEBP (lossless)
MIME Type                       : image/webp
Image Width                     : 736
Image Height                    : 515

Как и ожидалось, ничего полезного: ни EXIF, ни GPS, ни модели камеры.

b] Построение "отпечатка"

01_00

На изображении видно три массива суши:

P0: сам островок-курорт,
P1: остров справа,
P2: остров слева на переднем плане (с горной вершиной)

Построить корректную перспективную модель вида сверху не удалось — снимок явно сделан с дрона, и оценить высоту съёмки не получалось (в метаданных она тоже отсутствовала).

Поэтому расстояния и углы треугольника между тремя островами пришлось оценивать на глаз.

01_01

Для этого написан небольшой GUI 01_triangle_gui.py, который по клику записывает пиксельные координаты каждой точки по порядку и вычисляет геометрию треугольника.

Поскольку клик точно по центру объекта на глаз не идеален, к обоим значениям добавлен допуск ±20% при последующем поиске.

c] Поиск

Когда "отпечаток" зафиксирован, следующий шаг — сверить его с каждым реальным массивом суши на Земле.

В качестве набора данных использован раздельный набор земельных полигонов OpenStreetMap — land-polygons-split-4326, полный глобальный вектор береговой линии в WGS84 объёмом 882 МБ.

Эвристические фильтры (основанные исключительно на интуиции, без строгих доказательств) подбирались несколько дней методом проб и ошибок, пока не сложился рабочий набор правил.

01] Тропический пояс широт

$$ -30° \le latitude \le 30° $$

Остров на фото выглядит тропическим, поэтому всё, что находится вне тропического пояса, отбрасывается сразу, ещё до дорогостоящих геометрических вычислений.

Через этот фильтр проходят ровно 141 131 полигон суши.

02] Фильтр локальной плотности

$$ N_{5\text{km}}(p) \le 10 $$

$ N_{5\text{km}}(p) $ считает, сколько других центроидов попадает в радиус 5 км от точки (p). Порог — 10: если у островка больше 10 таких близких соседей, значит он расположен в плотном рифовом поле, на перегруженной береговой линии или в архипелаге с обилием мелких островов — а не в небольшой изолированной группе из 3-4 островов, как на фото.

Число кандидатов сократилось до 51 576.

03] Кластеризация

Для каждой оставшейся точки ищутся все остальные точки в радиусе 20 км (значение подобрано эвристически, на глаз по фотографии). Если у точки есть минимум 2 соседа на таком расстоянии (всего 3 точки), это кластер. Точки без кластера из 3+ соседей отбрасываются — они вообще не могут образовать треугольник.

tree = cKDTree(f_coords)
neigh = tree.query_ball_point(
                            f_coords, 
                            CLUSTER_RADIUS_KM / 111.0)
clusters = set(tuple(sorted(n)) for n in neigh if len(n) >= 3)

$$ \left|\{q : \text{dist}(p,q) \le 20\,\text{km}\}\right| \ge 3 $$

Итог — 23 500 кластеров.

04] Генерация троек

Для каждого кластера каждая комбинация из 3 точек внутри него становится треугольником-кандидатом. Это $ C(n, 3) $, число которых быстро растёт для больших кластеров: например, кластер из 60 точек сам по себе даёт 34 220 троек. Поэтому каждый кластер сначала ограничивается 60 точками, отобранными по размеру, а не случайно.

$$ \binom{n}{3} = \frac{n(n-1)(n-2)}{6} $$

def stratified_sample(idx_arr, area_arr, cap):
    order = np.argsort(area_arr[idx_arr])
    n_small = cap // 3
    n_large = cap // 3
    n_mid = cap - n_small - n_large
    mid_start = max(0, (len(idx_arr) - n_large - n_mid) // 2)
    keep = np.unique(np.concatenate([
        order[:n_small], 
        order[-n_large:], 
        order[mid_start:mid_start + n_mid],
    ]))
    return idx_arr[keep]

def gen_cluster_triples(idx_arr):
    local = np.array(list(
                itertools.combinations(range(len(idx_arr)), 3)), 
                dtype=np.int64)
    return idx_arr[local]

Выборка берёт треть маленьких островов, треть крупных и треть из середины распределения по размеру, вместо использования полного кластера или случайного среза.

23 500 кластеров дают в сумме 80 690 777 троек!

05] Сопоставление на GPU

Каждой тройке выделен один поток CUDA. Каждый поток сортирует свои 3 точки по площади суши, чтобы выделить P0 (наименьшую — островок-курорт), а затем по направлению обхода двух остальных точек определяет P1 и P2:

long long i = blockIdx.x * (long long)blockDim.x + threadIdx.x;
if (i >= n_triples) return;

int pos[3] = {0, 1, 2};
for (int a1 = 1; a1 < 3; a1++) 
{
    int key = pos[a1];
    double keyval = a[key];
    int j = a1 - 1;
    while (j >= 0 && a[pos[j]] > keyval) 
    {
        pos[j + 1] = pos[j];
        j--;
    }
    pos[j + 1] = key;
}

Выбор между P1 и P2 делается через векторное произведение в 2D, без ветвления по тому, из какого кластера пришла тройка — только по знаку:

$$ \text{cross} = x_a y_b - x_b y_a $$ $$ P1 = \begin{cases} a & \text{cross} > 0 \\ b & \text{cross} \le 0 \end{cases} $$

Проход от P0 к a, затем к b. Если cross > 0 — это поворот влево (против часовой стрелки), если cross < 0 — поворот вправо (по часовой). Тот же приём со знаком, что используется для определения направления изгиба трёх точек.

Далее вычисляются угол при P0 и отношение расстояний — по тем же формулам, что и на этапе построения "отпечатка", независимо каждым потоком:

$$ \theta_0 = \arccos\left(\frac{\vec{d_1} \cdot \vec{d_2}}{|\vec{d_1}||\vec{d_2}|}\right), \qquad r = \frac{|\vec{d_1}|}{|\vec{d_2}|} $$

Тройка проходит проверку, если угол, отношение расстояний, размер P0, расстояние между P0 и P1 и обе длины сторон укладываются в допуски "отпечатка". Прошедшие потоки записывают результат в общий выходной массив через атомарный счётчик, чтобы два потока, завершившиеся одновременно, не перезаписали друг друга:

if (hit) 
{
    unsigned long long slot = atomicAdd(out_count, 1ULL);
    out_p0[slot] = p0idx;
    out_p1[slot] = p1idx;
    out_p2[slot] = p2idx;
}

Результат работы ядра выводится прямо в консоль:

gpu: NVIDIA GeForce RTX 3050 (sm_86)
vram used: 5169 MB
kernel time: 204.1 ms

80,7 миллиона троек обрабатываются параллельно, по одному потоку на каждую. 158 784 проходят фильтр.

06] Дедупликация

Поскольку одна и та же физическая тройка может засветиться в нескольких потоках GPU, если она принадлежит более чем одному пересекающемуся кластеру, сырые совпадения сначала схлопываются по идентичности:

seen = set()
uniq = []
for i in range(len(p0_all)):
    key = (p0_all[i], p1_all[i], p2_all[i])
    if key not in seen:
        seen.add(key)
        uniq.append(i)

После дедупликации остаётся 8 915 уникальных троек.

07] Проверка на открытый прямоугольник

02_00

Каждая оставшаяся тройка проходит ещё одну проверку: действительно ли пространство рядом с ней — открытая вода, как показано на фото? Вдоль ребра P0→P1 строится прямоугольник, с той стороны, где не находится P2, и он проверяется на предмет наличия внутри чего-либо ещё из данных суши.

width = np.hypot(x1, y1)
u = np.array([x1, y1]) / width
v = np.array([-u[1], u[0]])

# p2 sits on the +v side by construction, 
# so the check goes on -v
length = 2 * width
corners_local = [
    (0, 0), (x1, y1),
    (x1 - v[0]*length, y1 - v[1]*length),
    (-v[0]*length, -v[1]*length),
]

Если что-либо, кроме трёх островов-кандидатов, пересекает этот прямоугольник, кандидат отбрасывается. Наличие там суши означает, что это не та открытая, ничем не загороженная вода, что видна на фото.

Из 8 915 уникальных троек остаётся 948.

Ниже — карта расположения этих 948 кандидатов.

02_01

d] Проверка формы кораллового острова

На этом этапе рассматривается только P0 — островок-курорт, и проверяется, действительно ли его форма похожа на коралловый остров.

1] `Компактность` — насколько форма близка к кругу:

Индекс Полсби — Поппера: $$ PP = \frac{4\pi \cdot \text{area}}{\text{perimeter}^2} $$

def compactness(row):
    return (4 * np.pi * row.area_km2) / (row.perim_km ** 2 + 1e-12)

03_00

1.0 соответствует идеальному кругу, чем ниже значение — тем более неровный или вытянутый контур. Коралловые острова обычно округлые за счёт волнового намыва, поэтому всё, что < 0.5, отбрасывается.

2] Проверка ореола из мелких островков

def micro_cay_count(gdf, sindex, lon, lat):
    dists_km = nearby.geometry.distance(pt) * 111.0
    mask = (dists_km > 0) 
           & (dists_km <= HALO_KM) 
           & (nearby["area_km2"].values < MICRO_KM2)
    return int(mask.sum())

Подсчитывались фрагменты суши площадью менее 0,05 км² в радиусе 1,5 км от P0 (эвристика). Настоящие рифовые системы разбрасывают вокруг основного острова мелкие песчаные отмели, а не представляют собой единственный изолированный массив суши (это стало понятно на собственном опыте). Поэтому требуется как минимум 1 такой фрагмент.

Обе проверки проходят 213 из 948 кандидатов.

e] Проверка овальной формы

Ещё один геометрический фильтр — по собственному полигону P0. Вокруг него строится минимальный повёрнутый прямоугольник, и по нему вычисляются два соотношения.

def aspect_and_fill(geom):
    mrr = geom.minimum_rotated_rectangle
    coords = list(mrr.exterior.coords)
    s1 = math.hypot(coords[1][0] - coords[0][0], 
                    coords[1][1] - coords[0][1])
    s2 = math.hypot(coords[2][0] - coords[1][0], 
                    coords[2][1] - coords[1][1])
    long_side, short_side = max(s1, s2), min(s1, s2)
    return long_side / short_side, geom.area / mrr.area

Соотношение сторон — длинная сторона к короткой стороне этого прямоугольника:

$$ \text{aspect} = \frac{\text{long side}}{\text{short side}} \in [1.05,\ 2.2] $$

Слишком близко к 1.0 — это практически идеальный круг, не совпадающий с чуть вытянутой формой на фото. Слишком высокое значение — форма вытянута сильнее, чем 2:1.

Коэффициент заполнения — какую часть ограничивающего прямоугольника занимает сама форма, и здесь есть точное тождество: любой эллипс заполняет ровно $ \pi / 4 $ площади своего минимального ограничивающего прямоугольника, независимо от степени вытянутости.

$$ \frac{\text{area}_{\text{ellipse}}}{\text{area}_{\text{box}}} = \frac{\pi}{4} \approx 0.785 $$

Это теоретический потолок для идеально гладкого овала. Реальные коралловые острова не являются идеальными эллипсами, поэтому порог отсечения задан как эвристическая безопасная доля от этого потолка:

$$ \text{FILL\_RATIO\_MIN} = 0.75 \times \frac{\pi}{4} \approx 0.589 $$

Форма должна сохранять как минимум 75% заполнения идеального эллипса, чтобы пройти проверку. Полумесяцы, кольца и изрезанные береговые линии оказываются значительно ниже этого порога, а сплошные округлые островки — нет.

Проверку проходят 137 из 213 кандидатов.

f] Проверка растительности по NDVI

Финальный этап связан с обращением к API, и он поставлен последним, поскольку упирается в скорость сети, а не вычислений.

Используется Earth Search от Element84 — публичный STAC API, индексирующий снимки Sentinel-2, размещённые в рамках программы открытых данных AWS, бесплатно и без ключа API.

Ознакомиться с ним можно по адресу https://earth-search.aws.element84.com/v1

Проверяется, действительно ли P0 покрыт растительностью — пальмами, а не голым песком или скалой. Из публичного каталога STAC берётся самый свежий снимок Sentinel-2 с низкой облачностью над этой точкой, и по нему считываются значения красного и ближнего инфракрасного каналов в точной пиксельной позиции.

$$ \text{NDVI} = \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}} $$

Живая растительность сильно отражает в ближнем инфракрасном диапазоне и поглощает красный свет, поэтому здоровый пальмовый покров даёт NDVI значительно выше 0, а голый песок или открытая вода — около 0 или ниже.

04_00

Иллюстрация взята из блога Geoawesome.

Порог установлен на уровне 0.6 — достаточно высоко, чтобы требовать реального древесного покрова, а не отдельных разрозненных объектов.

Проверку по NDVI проходят 66 из 137.

g] Проверка рельефа и наличия гор

05_00

Последняя проверка перед итоговым результатом. Здесь два условия:

сам P0 должен быть низким и плоским, что соответствует небольшому рифовому островку,
у P2 должен быть реальный возвышенный рельеф в том направлении, куда была направлена камера.

"Фронт" снимка — биссектриса между азимутом на P1 и азимутом на P2:

$$ \theta(P_0, P_i) = $$ $$ \text{atan2}\Big(\sin(\Delta\lambda)\cos\phi_i,\ \cos\phi_0\sin\phi_i - \sin\phi_0\cos\phi_i\cos(\Delta\lambda)\Big) $$

$$ \theta_{\text{front}} = $$ $$ \theta(P_0, P_2) + \frac{\big((\theta(P_0,P_1) - \theta(P_0,P_2) + 180) \bmod 360\big) - 180}{2} $$

Это даёт единый азимут — направление, куда был направлен объектив. От него веером разворачиваются точки выборки в диапазоне ±50° вокруг этого азимута, на удалении от 2 до 20 км:

$$ (\text{lat}, \text{lon}) = \Big(\text{lat}_0 + \frac{r\cos\theta}{111},\ \ \text{lon}_0 + \frac{r\sin\theta}{111\cos(\text{lat}_0)}\Big) $$

Каждая из этих точек сверяется с реальными тайлами 30-метрового Copernicus DEM.

Copernicus DEM GLO-30, публикуемый программой Copernicus Евросоюза, размещён в виде бесплатных облачно-оптимизированных GeoTIFF в рамках AWS Open Data, без учётной записи и ключа.

Подробнее — https://registry.opendata.aws/copernicus-dem/

В итоге выживание кандидата определяют два простых эвристических условия (да, к этому моменту всё уже стало эвристикой):

$$ \text{elev}(P_0) \le 50\text{m} $$ $$ 100\text{m} \le \max_{\text{arc}}(\text{elev}) \le 500\text{m} $$

05_01

На этой схематичной иллюстрации пунктирная линия — азимут "фронта" камеры, а сектор — веер поиска ±50°, развёрнутый до 20 км для проверки рельефа.

Проверку по рельефу проходят 26 из 66.

Все 26 оставшихся кандидатов расположены в Южной Азии, Австралии и Океании, кроме одного — рядом с Бразилией!

05_02

h] Итоговый отчёт

На последнем этапе финальные кандидаты приводятся в вид, удобный для проверки глазами. Каждому выжившему кандидату присваивается название страны через точечный поиск в полигоне по файлу государственных границ, а также прямая ссылка на спутниковый снимок Google Maps для P0, P1 и P2.

Результат — обычная HTML-таблица: индекс, страна, три кликабельные пары координат в строке.

06_00

По этому итоговому списку каждый кандидат проверялся визуально.

Первые 7 строк оказались явно не подходящими.

06_01

Но восьмая строка в таблице привела в страну Микронезия (впервые пришлось узнать, что существует страна с таким названием):

06_03

и подтверждение пришло через P1 и P2:

06_04

Это и есть решение...

Посмотреть на карте можно здесь: Google Maps

i] Итоговые ответы

a) What is the name of the resort?

$$ \text{Oan} $$

b) What are the coordinates of the island?

$$7^\circ\,21^\prime\,48.4^{\prime\prime}\,\text{N} \qquad 151^\circ\,45^\prime\,20.7^{\prime\prime}\,\text{E}$$

$$ \text{или} $$

$$7.363444^\circ,\ 151.755750^\circ$$

c) In which cardinal direction was the 
camera facing when the photo was taken?

$$ \because\quad \theta = \text{atan2}\Big(\sin(\Delta\lambda)\cos\phi_1,\ \cos\phi_0\sin\phi_1 - \sin\phi_0\cos\phi_1\cos(\Delta\lambda)\Big) $$

$$ P_0 = (7.3633,\ 151.755983), \quad P_1 = (7.386573,\ 151.739534) $$

$$ \therefore\quad \theta = 324.97^\circ \implies \textbf{СЗ} $$

j] Данные и лицензии

Полигоны береговой линии:
land-polygons-split-4326 © участники OpenStreetMap, доступно по лицензии Open Database License (ODbL) 1.0. Наборы кандидатов и итоговый отчёт в репозитории являются производной базой данных и публикуются под той же лицензией.

Спутниковые снимки:
Содержит модифицированные данные Copernicus Sentinel 2025-2026, полученные через Earth Search компании Element 84 в рамках AWS Open Data.

Границы стран:
Natural Earth 10m admin-0, общественное достояние.

Задача и исходное фото:
OSINT Exercise #004 от Sofia Santos (gralhix).

Спутниковые скриншоты в разделе (h) взяты из Google Maps / Google Earth.