Как найти остров по одной фотографии: геолокация через геометрию и CUDA

Перед вами фотография островного курорта.

1. Как называется курорт?

2. Какие координаты у острова?

3. В какую сторону света была направлена камера в момент съёмки?

На мой взгляд, пытаться ответить на эти вопросы через google lens – значит упустить возможность повеселиться. Поэтому я взялся за задачу с помощью математики и программирования.

Метаданные

Разумеется, первым делом смотрим метаданные. Запустил это на своём Void Linux:

> exiftool main.png

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

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

Строим «отпечаток»

На фото видно три участка суши:

• P0: сам островок;

• P1: остров справа;

• P2: остров слева на переднем плане (с горной вершиной).

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

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

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

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

Поиск

Отпечаток зафиксирован – теперь нужно сверить с ним каждый реальный участок суши на Земле!

В качестве датасета я взял разбитый на части набор полигонов суши из OpenStreetMap – land-polygons-split-4326 ^[1]. Это полные векторы береговых линий всего мира в WGS84 размером 882 МБ.

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

Ограничивающая рамка по тропическим широтам

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

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

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

считает, сколько других центроидов попадает в радиус 5 км от точки (p). Предел – 10: если у островка больше 10 соседей на таком расстоянии, значит, он сидит в плотном рифовом поле, у изрезанного побережья или в скоплении архипелага – а не в маленькой изолированной группе из 3-4 островов, как на фото.

Это сократило число кандидатов до 51 576.

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

Для каждой уцелевшей точки ищем все остальные точки в радиусе 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)

В итоге остаётся 23 500 кластеров.

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

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

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 троек!!

Сопоставление на 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 определяются через двумерное векторное произведение – без всякого ветвления по тому, из какого кластера пришла тройка, только по знаку:

Идём от P0 к a, затем к b. Если cross > 0 – это поворот налево (против часовой стрелки). Если cross < 0 – поворот направо (по часовой стрелке). Это тот же приём со знаком, которым определяют, в какую сторону изгибаются три точки.

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

Тройка проходит, если угол, отношение, размер 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 млн троек, каждая – в своём потоке, параллельно. Маску проходят 158 784.

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

Одна и та же физическая тройка может попасть сразу в несколько 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 уникальных троек.

Открытый прямоугольник

Каждая уцелевшая тройка проходит ещё одну проверку: действительно ли рядом открытая вода, как на фото? Вдоль ребра 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),
]

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

8 915 уникальных троек сократились до 948.

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

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

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

Компактность – насколько форма близка к кругу:

Индекс Полсби-Поппера (Polsby-Popper Score):

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

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

Проверка ореола из микро-островков:

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 кандидатов.

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

Ещё один геометрический фильтр по самому полигону 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 ratio) – это длинная сторона рамки, делённая на короткую:

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

Коэффициент заполнения (fill ratio) показывает, какую часть ограничивающей рамки реально занимает форма, и за ним стоит точное тождество: любой эллипс заполняет ровно своего минимального ограничивающего прямоугольника, как бы его ни растягивали.

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

FILL_RATIO_MIN

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

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

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

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

Мы подключимся к Earth Search от Element84 – это публичное STAC API, которое индексирует снимки Sentinel-2, размещённые в программе AWS Open Data. Бесплатно и без API-ключа.

Посмотреть можно здесь: https://earth-search.aws.element84.com/v1 ^[2]

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

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

Это изображение я взял из блога Geoawesome. ^[3]

Порог выставлен на 0.6 – достаточно высоко, чтобы требовать реального древесного покрова, а не отдельных пятен.

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

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

Последняя проверка перед финалом. Условий два:

• сам P0 должен быть низким и плоским, как и подобает маленькому рифовому островку,

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

«Перед» кадра – это биссектриса между азимутом на P1 и азимутом на P2:

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

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

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

Подробнее – здесь: https://registry.opendata.aws/copernicus-dem/ ^[4]

И наконец, судьбу кандидата решают два простых эвристических условия:

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

Проверку по высоте проходят 26 из 66.

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

Финальный отчёт

И вот, последний этап – он просто делает финальных кандидатов пригодными для проверки глазами. Для каждого уцелевшего кандидата смотрим, в какую страну попадает точка (проверкой «точка в полигоне» по файлу госграниц), а затем формируем прямую ссылку на спутниковые снимки Google Maps для P0, P1 и P2.

На выходе – обычная HTML-таблица: индекс, страна и три кликабельные пары координат в каждой строке.

Получился такой финальный список – давайте проверим каждый пункт глазами. Разбирать по одному я тут не буду, но первые 7 для меня совсем мимо.

…пока я не открыл восьмой пункт – из страны Микронезия (впервые узнал, что есть страна с таким названием):

и убедился в этом по P1 и P2:

Вот оно, решение! Посмотреть можно здесь, на google maps ^[5]

Наконец, ответы!

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

2) Какие координаты у острова?

3) В какую сторону света была направлена камера в момент съёмки?

Посмотреть, склонировать и запустить у себя все файлы с кодом, а также финальный отчёт со всеми инструкциями можно здесь, на GitHub. ^[6]