Таинственный остров: находим геолокацию с помощью геометрии и программирования GPU CUDA / Хабр

Свой пост я написал после участия в соревнованиях gralhix 004, организованных Софией Сантос | Gralhix.

Задача

Это фотография островного курорта.

Вопросы:

а) Как называется курорт?
б) Каковы координаты острова?
в) В какую сторону света была направлена камера, когда делали снимок?

На мой взгляд, решение этой задачи при помощи Google Объектива будет потраченной впустую возможности развлечься, поэтому я захотел решить её при помощи математики и программирования.

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

Разумеется, первым делом следует посмотреть метаданные. Я выполнил в Void Linux следующую команду:

> exiftool main.png

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

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

б] Собираем фингерпринт

На фотографии есть три массива суши:

P0: сам островок,
P1: остров справа,
P2: остров слева (с горой)

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

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

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

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

в] Поиск

Определившись с фингерпринтом, можно приступать к сверке с ним всех массивов суши на Земле!

Я использовал в качестве датасета split land polygon set OpenStreetMap (land‑polygons‑split-4326), в нём хранятся все векторные береговые линии Земли в формате WGS84; весь датасет занимает 882 МБ.

Я создал эвристические фильтры (интуитивно, без математических доказательств), потратив несколько дней на настройку значений путём проб и ошибок, пока не получил готовые фильтры.

01] Ограничивающий прямоугольник тропической широты

−30° ≤ широта ≤ 30°

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

Через этот полосовой фильтр прошёл ровно 141131 полигон суши.

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

$N_{5km}(p)≤10$

$N_{5km}(p)≤10$ определяет, сколько центров масс находится в пределах 5 километров от точки (p). Максимальное значение равно 10: если у острова больше десяти соседей на таком расстоянии, то он находится на плотном рифе, заполненной островами береговой линии или в архипелаге, а не в маленькой группе из 3–4 островов, как мы видим на фото.

Этот фильтр уменьшил количество кандидатов до 51576.

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

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

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$

Так мы снизили количество до 23500 кластеров.

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

В каждом кластере все комбинации трёх входящих в него точек становятся треугольником‑кандидатом. Это, что приводит к комбинаторному взрыву для больших кластеров; например, кластер из 60 точек уже даёт 34220 троек. Поэтому каждый кластер сначала ограничивается 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]

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

23500 суммарно дают 80690777 троек!

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

Я выделил каждой тройке поток CUDA. Каждый поток сортирует свои три точки по площади суши для выбора 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 происходит на основе двухмерного векторного произведения: никакого ветвления в зависимости от того, из какого кластера взята тройка, берётся только знак:

$\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;
}

Вывод в CLI напрямую из ядра:

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

80,7 миллиона троек на входе, по одной на поток. Через маску прошли 158784.

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)

После дедупликации осталось 8915 уникальных троек.

07] Свободный прямоугольник

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

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

# p2 по построению находится на стороне +v, 
# поэтому проверка выполняется для -v
length = 2 * width
corners_local = [
    (0, 0), (x1, y1),
    (x1 - v[0]*length, y1 - v[1]*length),
    (-v[0]*length, -v[1]*length),
]

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

Из 8915 уникальных троек осталось 948.

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

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

На этом этапе мы рассматриваем только 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)

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 (это просто моя эвристика). Реальные рифовые системы — это не одна изолированная масса суши, они разбрасывают вокруг основного острова крошечные песчаные косы. Поэтому нам нужна хотя бы одна.

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

д] Проверка овальности

Ещё один геометрический фильтр полигона 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

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

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

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

Соотношение Fill — показывает, какую долю ограничивающего прямоугольника заполняет форма острова: любой эллипс заполняет ровно π/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 кандидатов.

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

Мы дошли до последнего этапа API: я поместил его в конец, потому что его основным ограничением будут не вычислительные ресурсы, а сетевые.

Мы подключимся к публичному STAC API Earth Search компании Element84, индексирующему спутниковые снимки Sentinel-2, которые хостятся в программе AWS Open Data бесплатно и без ключа 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 или отрицателен.

Можете посмотреть это изображение, которое я взял из блога Geoawesome.

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

Проверку NDVI прошло 66 из 137 кандидатов.

ж] Проверка высоты и горы

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

Сама 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)$

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

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

Подробную информацию можно изучить здесь: 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}$

На этом абстрактном чертеже пунктирная линия — это направление камеры, сектор — это дуга поиска ±50° с радиусом 20 км для поиска возвышений.

Проверку на высоту прошли 26 из 66 кандидатов.

Ниже показаны все оставшиеся кандидаты, расположенные в Южной Азии, Австралии и Океании. Только одна из точек находится рядом с Бразилией.

з] Окончательный отчёт

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

Вывод представляет собой HTML‑таблицу, в каждой строке которой указаны индекс, страна и три ссылки на координаты.

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

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

А потом я открыл восьмой пункт таблицы, находящийся в Микронезии (так я узнал о стране под названием Микронезия):

Затем проверил P1 и P2:

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

Можете сами посмотреть его на Google Картах

и] И, наконец, ответы на вопросы...

a) Как называется курорт?

Oan

б) Каковы координаты острова?

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

или

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

в) В какую сторону света была направлена камера, когда делали снимок ?

$\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{NW}$

Северо‑запад

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

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

Высоты:
Copernicus DEM GLO-30. © DLR e.V. 2010–2014 и © Airbus Defence and Space GmbH 2014–2018, предоставленные в рамках программы COPERNICUS Евросоюзом и ESA; с сохранением всех прав.

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

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

Challenge & source photo :
OSINT Exercise #004, София Сантос (gralhix).

Скриншоты спутниковых фотографий из раздела (h) взяты из Google Карт / Google Планеты Земля

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

Таинственный остров: находим геолокацию с помощью геометрии и программирования GPU CUDA