Geolocating a random island using geometry and CUDA programming

Aug 19, 2026 07:19 PM - 2 hours ago 4

16-08-2026

NOTE: this is simply a genuine quality work, didnt usage LLM generation.

I'm penning this page arsenic a writeup for this situation gralhix 004 made by Sofia Santos | Gralhix.

You tin view, clone and locally effort each codification files and the last study pinch each instructions here astatine github.

Task briefing:

main

This is simply a photograph of a edifice located connected an island.

a) What is the sanction of the resort?
b) What are the coordinates of the island?
c) In which cardinal guidance was the camera facing erstwhile the photograph was taken?

In my opinion, solving this situation pinch google lens is wasting a nosy opportunity, truthful decided to lick it pinch mathematics and programming.

a] Metadata

Of course, first point u look for is the metadata. Ran that connected my linux void:

> exiftool main.png File Type : WEBP (lossless) MIME Type : image/webp Image Width : 736 Image Height : 515

As expected, thing useful here. No EXIF, nary GPS, nary camera make aliases model.

b] Building the fingerprint

01_00

U tin spot from the img, location are 3 landmasses:

P0: the islet itself,
P1: the correct island,
P2: the near beforehand land ( having upland highest )

I couldnt make a correct position exemplary of birdview of this image, arsenic intelligibly the image is taken by a drone and cant estimate the elevation astatine each (and not recovered successful the metadata).

So I had to estimate that by intuition, I conscionable want the comparative distances betwixt the 3 islands and angles of that triangle.

01_01

I built a mini click GUI 01_triangle_gui.py that records pixel coordinates for each constituent successful bid and computes the triangle's geometry.

Since clicking nonstop centers by oculus isn't perfectly precise, I added a ±20% tolerance set astir some values erstwhile searching.

c] SEARCH

With the fingerprint locked in, the adjacent measurement is checking each existent landmass connected Earth against it !

I utilized OpenStreetMap's divided onshore polygon group arsenic the dataset land-polygons-split-4326, afloat world coastline vectors successful WGS84 which has size of 882 MB.

I created heuristic filters (all by conscionable intuition and non tangible proofs), spent days (yea afloat days) tweaking values and tons of proceedings and correction 😭 untill I sewage this moving filters recipe.

01] Tropical latitude bounding box

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

the islet successful the photograph sounds arsenic tropical, truthful I decided that thing extracurricular the tropics is thrown retired immediately, earlier doing immoderate costly geometry work.

Exactly 141,131 onshore polygons past that set filter.

02] Local density filter

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

$ N_{5\text{km}}(p) $ counts really galore different centroids autumn wrong 5km of constituent (p). Cap is 10: if an islet has much than 10 neighbors that close, it's sitting successful a dense reef field, a crowded coastline aliases a archipelago clutter, not a mini isolated 3-4 land group for illustration the photograph shows.

This dropped candidates down to 51,576.

03] Clustering

For each surviving point, find each different constituent wrong 20km (heuristic, by oculus from the image). If it has astatine slightest 2 neighbors that adjacent (3 points total), it's a cluster. Points pinch nary cluster of 3+ adjacent are dropped, they can't shape a triangle astatine all.

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 $$

That collapses down to 23,500 clusters.

04] Generating Triplets

For each cluster, each operation of 3 points wrong it becomes a campaigner triangle. That's $ C(n, 3) $, which explodes accelerated for large clusters, for example: a cluster of 60 points already gives 34,220 triples connected its own. So each cluster gets capped astatine 60 points first, sampled by size, not randomly.

$$ \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]

The sampling takes a 3rd mini islands, a 3rd large, a 3rd from the mediate of the size distribution, alternatively of the afloat cluster aliases a random cut.

23,500 clusters nutrient 80,690,777 triples full !!

05] Matching, connected the GPU

I gave each triple 1 CUDA thread. Each thread sorts its 3 points by onshore area to prime retired P0 (smallest, the edifice islet), past uses the winding guidance of the different 2 to delegate P1 and 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 vs P2 comes from a 2D transverse product, nary branching connected which cluster the triple came from, conscionable the sign:

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

Walk from P0 to a, past to b. If transverse > 0, that's a near move (counterclockwise). If transverse < 0, it's a correct move (clockwise). It's the aforesaid motion instrumentality utilized to show if 3 points curve 1 measurement aliases the other.

then perspective astatine P0 and the region ratio, aforesaid formulas arsenic the fingerprint step, computed independently by each thread:

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

A triple survives if angle, ratio, P0's size, the separation betwixt P0 and P1, and some broadside lengths each onshore wrong the fingerprint's tolerance windows. Threads that walk constitute their consequence into a shared output array utilizing an atomic counter, truthful 2 threads finishing astatine the aforesaid clip ne'er overwrite each other:

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

Now printed successful the CLI straight from the kernel:

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

80.7 cardinal triples spell in, 1 thread each, successful parallel. 158,784 walk the mask.

06] Dedup

Since aforesaid beingness triple tin get deed by aggregate GPU threads if it belonged to much than 1 overlapping cluster, truthful earthy matches get collapsed by personality first:

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 unsocial triples aft dedup.

07] The Open Rectangle

02_00

Every surviving triple gets 1 much test: is the abstraction adjacent to it really unfastened water, for illustration the photograph shows ? A rectangle gets built on the P0→P1 edge, connected whichever broadside P2 is not on, past checked against the onshore dataset for thing other sitting wrong it.

width = np.hypot(x1, y1) u = np.array([x1, y1]) / width v = np.array([-u[1], u[0]]) # p2 sits connected the +v broadside by construction, # truthful the cheque goes connected -v length = 2 * width corners_local = [ (0, 0), (x1, y1), (x1 - v[0]*length, y1 - v[1]*length), (-v[0]*length, -v[1]*length), ]

If thing different than the 3 campaigner islands themselves intersects that rectangle, the campaigner is dropped. Land sitting location intends it's not the open, unobstructed h2o the photograph really shows.

8,915 unsocial triples down to 948.

and beneath is the representation of places of the 948 candidates.

02_01

d] Coral Cay Shape Check

In this stage, we look only astatine P0, the edifice islet, and cheque whether its style really looks for illustration a coral cay.

1] Compactness, really adjacent to a circle the style is:

Polsby Popper Score: $$ 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 is simply a cleanable circle, little intends a much jagged aliases elongated outline. Coral cays thin to beryllium information from activity deposition, truthful thing < 0.5 gets dropped.

2] Micro Cay Halo Check:

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())

We Count onshore fragments nether 0.05 km² wrong 1.5km of P0 ( conscionable heuristic ). Real reef systems scatter mini sandbars astir the main island, not conscionable 1 isolated landmass (I knew that pinch the hardway 😭). So we request astatine slightest 1.

213/948 candidates past some checks.

e] Oval Shape Check

Another geometric select connected P0's ain polygon. Fits the minimum rotated rectangle astir it and measures 2 ratios from that box.

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 is agelong broadside complete short broadside of that box:

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

Too adjacent to 1.0 and it's fundamentally a cleanable circle, not the somewhat elongated style successful the photo. Too precocious are shapes excessively overmuch elongated much than 2:1.

Fill ratio is really overmuch of that bounding container the style really fills, and this 1 has an personality down it: immoderate ellipse fills precisely $ \pi / 4 $ of its ain minimum area bounding rectangle, sloppy of really stretched it is.

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

that's the theoretical ceiling for a perfectly soft oval. Real coral cays aren't cleanable ellipses, truthful the cutoff is group arsenic a heuristic safe fraction of that ceiling:

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

A style needs to clasp astatine slightest 75% of a cleanable ellipse's capable to survive. Crescents, rings, and notched coastlines autumn good beneath that, coagulated rounded cays don't.

137/213 candidates survive.

f] NDVI Vegetation Check

We reached the last API phase, I put it astatine the end, because it is web bound not compute bound.

We gonna link to Earth Search, tally by Element84, a nationalist STAC API that indexes Sentinel-2 imagery hosted connected AWS's Open Data program, free, nary API key.

You tin look astatine it https://earth-search.aws.element84.com/v1

We now cheque whether P0 is really vegetated, thenar cover, not bare soil aliases rock. It pulls the astir caller debased unreality Sentinel-2 segment complete the constituent from a nationalist STAC catalog, samples the reddish and adjacent infrared bands astatine that nonstop pixel.

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

Live vegetation reflects powerfully successful adjacent infrared and absorbs reddish light, truthful patient thenar screen pushes NDVI good supra 0, bare soil aliases unfastened h2o sits adjacent 0 aliases negative.

04_00

You tin position this image I sewage from this bully Geoawesome Blog.

Threshold is group astatine 0.6, precocious capable to require existent character cover, not conscionable scattered units.

66/137 past the NDVI check.

g] Elevation & Mountain Check

05_00

Last cheque earlier the last reveal. There are 2 conditions:

P0 itself must beryllium debased and flat, accordant pinch a mini reef islet,
P2 must person existent elevated terrain successful the guidance the camera was really facing.

The "front" of the changeable is the bisector betwixt the base to P1 and the base to 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} $$

That gives 1 heading, the guidance the lens was pointed. From there, a instrumentality of sample points gets swept ±50° astir that heading, astatine radii from 2km retired to 20km:

$$ (\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) $$

Every 1 of those points gets sampled against existent 30m Copernicus DEM tiles.

Copernicus DEM GLO-30, published by the EU's Copernicus program, hosted arsenic free nationalist Cloud-Optimized GeoTIFFs connected AWS Open Data, nary relationship aliases cardinal needed.

For much info, you tin position https://registry.opendata.aws/copernicus-dem/

Finally, those 2 elemental heuristic conditions determine endurance (yea I know, everything became heuristic haha):

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

05_01

You tin spot from this absurd chart image, the dashed statement is the camera's beforehand bearing, the wedge is the ±50° hunt arc swept retired to 20km for the elevation check.

26/66 past the elevation check.

You tin spot the 26 survivors, each are located successful confederate Asia, Australia and Oceania, isolated from 1 adjacent Brazil!

05_02

h] Final Report

Finally, past stage, it conscionable makes the last candidates checkable by eye. Each subsister gets its state sanction via a constituent successful polygon lookup against a state bound file, past a nonstop Google Maps outer nexus for P0, P1, and P2.

Output is simply a plain HTML table, index, country, 3 clickable coordinate pairs per row.

06_00

I sewage this last list, lets cheque each 1 by eye.

Won't spell 1 by 1 here, but those first 7 are wholly disconnected for me.

06_01

Till I opened that 8th 1 successful the array of state of Micronesia 😍 (first clip to cognize that a state named Micronesia):

06_03

and ensured done P1 and P2:

06_04

and that is the solution 🥳 ...

you tin position it present connected google maps

i] FINALLY, ANSWERS ...

a) What is the sanction 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{or} $$

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

c) In which cardinal guidance was the camera facing erstwhile the photograph 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{NW} $$

j] Data & Licenses

Coastline polygons:
land-polygons-split-4326 © OpenStreetMap contributors, disposable nether the Open Database License (ODbL) 1.0. The campaigner sets and last study successful the repo are a Derived Database and are published nether the aforesaid license.

Elevation :
Copernicus DEM GLO-30. © DLR e.V. 2010-2014 and © Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by the European Union and ESA; each authorities reserved.

Satellite imagery :
Contains modified Copernicus Sentinel information 2025-2026, accessed done Earth Search by Element 84 connected AWS Open Data.

Country boundaries :
Natural Earth 10m admin-0, nationalist domain.

Challenge & root photograph :
OSINT Exercise #004 by Sofia Santos (gralhix).

Satellite screenshots successful conception (h) are from Google Maps / Google Earth