forked from
anil.recoil.org/monopam-myspace
My aggregated monorepo of OCaml code, automaintained
Fix GeoTessera overlay alignment by reprojecting tiles from UTM to WGS84
The GeoTessera embedding tiles have pixel grids in UTM projection, but the
notebook was overlaying them as if they were on a regular lat/lng grid. This
caused ~3-6% east-west distortion at mid-latitudes (e.g. ~440m at 52N).
- Add pure-OCaml UTM projection module (utm.ml) with WGS84 ellipsoid
forward/inverse transforms
- Add reproject_tile function that resamples each tile from its native UTM
grid onto a regular WGS84 grid using nearest-neighbor interpolation
- Mosaic now reprojects all tiles before assembly, matching the Python
reference implementation (ucam-eo/tessera-interactive-map) which uses
rasterio.warp.reproject for the same purpose
- Return computed WGS84 bounds from fetch_mosaic_sync instead of requiring
callers to hardcode snap ± 0.05
- Fix jtw relativize_or_fallback for dune exec paths (dynamic_cmis.json
generation was silently skipped for locally-built packages)
- Consolidate deploy-site.sh as thin wrapper around build-site.sh
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
+274
-189
+6
-5
build-site.sh
+6
-5
build-site.sh
···
40
40
echo "=== Step 1: Build and register plugins ==="
41
41
cd "$MONO"
42
42
dune build @install
43
-
dune install 2>/dev/null
43
+
dune install --force 2>&1 | grep -v "Installing\|Overwriting" || true
44
44
echo " plugins registered"
45
45
46
46
echo ""
···
59
59
cp -rf "$DUNE_DOC/reference/"* "$SITE/reference/"
60
60
mkdir -p "$SITE/odoc.support/"
61
61
cp -rf "$DUNE_DOC/odoc.support/"* "$SITE/odoc.support/"
62
+
cp -f "$DUNE_DOC/sidebar.json" "$SITE/sidebar.json" 2>/dev/null || true
62
63
echo " assembled into $SITE/"
63
64
64
65
echo ""
···
66
67
if [ -f "$SITE/_opam/worker.js" ]; then
67
68
echo " universe already exists, skipping (use --fresh to rebuild)"
68
69
else
69
-
jtw opam astring base brr note mime_printer fpath rresult \
70
+
dune exec -- jtw opam astring base brr note mime_printer fpath rresult \
70
71
opam-format bos odoc.model tyxml yojson uri jsonm \
71
72
js_top_worker-widget-leaflet \
72
73
tessera-geotessera-jsoo tessera-viz-jsoo \
···
97
98
trap 'rm -rf "$UNIVERSES"' EXIT
98
99
99
100
echo " building default universe (cmdliner, 5.2.0+ox switch)..."
100
-
jtw opam --switch=5.2.0+ox -o "$UNIVERSES/default" cmdliner
101
+
dune exec -- jtw opam --switch=5.2.0+ox -o "$UNIVERSES/default" cmdliner
101
102
102
103
echo " building v3 universe (cmdliner, 5.2.0+ox switch)..."
103
-
jtw opam --switch=5.2.0+ox -o "$UNIVERSES/v3" cmdliner
104
+
dune exec -- jtw opam --switch=5.2.0+ox -o "$UNIVERSES/v3" cmdliner
104
105
105
106
echo " building oxcaml universe (5.2.0+ox switch)..."
106
-
jtw opam --switch=5.2.0+ox -o "$UNIVERSES/oxcaml"
107
+
dune exec -- jtw opam --switch=5.2.0+ox -o "$UNIVERSES/oxcaml"
107
108
108
109
for d in universe universe-v2 universe-v3 universe-oxcaml; do
109
110
rm -rf "$DEMO_DIR/$d"
+10
-132
deploy-site.sh
+10
-132
deploy-site.sh
···
1
1
#!/bin/bash
2
-
# Build and optionally serve the full jon.recoil.org site.
2
+
# Build and serve the full jon.recoil.org site.
3
3
#
4
-
# Dune outputs go into _build/ (which dune controls and may wipe).
5
-
# We assemble the final site into _site/ so that expensive artifacts
6
-
# like universes persist across rebuilds.
4
+
# This is a convenience wrapper around build-site.sh that defaults to
5
+
# serving on port 8080 after building.
7
6
#
8
7
# Usage:
9
8
# ./deploy-site.sh # build everything and serve on port 8080
···
13
12
set -euo pipefail
14
13
15
14
MONO=$(cd "$(dirname "$0")" && pwd)
16
-
SITE="$MONO/_site"
17
-
DUNE_SITE="$MONO/_build/default/site/_html"
18
-
DUNE_DOC="$MONO/_build/default/_doc/_html"
15
+
16
+
# Translate --no-serve → omit --serve; otherwise pass --serve plus any other flags.
17
+
BUILD_ARGS=()
19
18
SERVE=true
20
-
FRESH=false
21
-
22
19
for arg in "$@"; do
23
20
case "$arg" in
24
21
--no-serve) SERVE=false ;;
25
-
--fresh) FRESH=true ;;
22
+
*) BUILD_ARGS+=("$arg") ;;
26
23
esac
27
24
done
28
25
29
-
if $FRESH; then
30
-
echo "=== --fresh: removing _site ==="
31
-
rm -rf "$SITE"
32
-
fi
33
-
34
-
mkdir -p "$SITE"
35
-
36
-
# Ensure we're on the right switch.
37
-
export OPAMSWITCH=5.2.0+ox
38
-
eval "$(opam env)"
39
-
40
-
echo "=== Step 1: Build and register plugins ==="
41
-
cd "$MONO"
42
-
dune build @install
43
-
echo " plugins registered"
44
-
45
-
echo ""
46
-
echo "=== Step 2: Build site content ==="
47
-
dune build @site 2>&1 | grep -v '^Warning\|^File\|^$' | tail -5 || true
48
-
echo " dune @site done"
49
-
50
-
echo ""
51
-
echo "=== Step 3: Build reference docs ==="
52
-
dune build @doc 2>&1 | tail -5 || true
53
-
echo " dune @doc done"
54
-
55
-
echo ""
56
-
echo "=== Step 4: Assemble site ==="
57
-
# Copy site HTML from dune output into _site.
58
-
# rsync would be ideal but cp -rf works fine.
59
-
cp -rf "$DUNE_SITE/"* "$SITE/"
60
-
# Merge reference docs on top (site has reference/index.html, @doc adds packages).
61
-
cp -rf "$DUNE_DOC/reference/"* "$SITE/reference/"
62
-
mkdir -p "$SITE/odoc.support/"
63
-
cp -rf "$DUNE_DOC/odoc.support/"* "$SITE/odoc.support/"
64
-
# sidebar.json is fetched at runtime by the sidebar JS
65
-
cp -f "$DUNE_DOC/sidebar.json" "$SITE/sidebar.json"
66
-
echo " assembled into $SITE/"
67
-
68
-
echo ""
69
-
echo "=== Step 5: Build site universe (notebooks) ==="
70
-
if [ -f "$SITE/_opam/worker.js" ]; then
71
-
echo " universe already exists, skipping (use --fresh to rebuild)"
72
-
else
73
-
jtw opam astring base brr note mime_printer fpath rresult \
74
-
opam-format bos odoc.model tyxml yojson uri jsonm \
75
-
js_top_worker-widget-leaflet \
76
-
tessera-geotessera-jsoo tessera-viz-jsoo \
77
-
onnxrt -o "$SITE/_opam"
78
-
echo " universe built → $SITE/_opam/"
79
-
fi
80
-
81
-
echo ""
82
-
echo "=== Step 6: Deploy onnxrt example assets ==="
83
-
SENTIMENT_SRC="$MONO/onnxrt/example/sentiment"
84
-
if [ ! -f "$SENTIMENT_SRC/model_quantized.onnx" ]; then
85
-
echo " downloading DistilBERT model..."
86
-
bash "$SENTIMENT_SRC/download_model.sh"
26
+
if $SERVE; then
27
+
BUILD_ARGS+=("--serve")
87
28
fi
88
-
cp "$SENTIMENT_SRC/vocab.txt" "$SITE/_opam/vocab.txt"
89
-
cp "$SENTIMENT_SRC/model_quantized.onnx" "$SITE/_opam/model_quantized.onnx"
90
-
cp "$MONO/onnxrt/example/add.onnx" "$SITE/_opam/add.onnx"
91
-
echo " deployed vocab.txt + model_quantized.onnx + add.onnx → $SITE/_opam/"
92
29
93
-
echo ""
94
-
echo "=== Step 7: Build demo universes ==="
95
-
DEMO_DIR="$SITE/reference/odoc-interactive-extension"
96
-
97
-
if [ -f "$DEMO_DIR/universe/worker.js" ]; then
98
-
echo " demo universes already exist, skipping (use --fresh to rebuild)"
99
-
else
100
-
UNIVERSES=$(mktemp -d)
101
-
trap 'rm -rf "$UNIVERSES"' EXIT
102
-
103
-
echo " building default universe (cmdliner, 5.2.0+ox switch)..."
104
-
jtw opam --switch=5.2.0+ox -o "$UNIVERSES/default" cmdliner
105
-
106
-
echo " building v3 universe (cmdliner, 5.2.0+ox switch)..."
107
-
jtw opam --switch=5.2.0+ox -o "$UNIVERSES/v3" cmdliner
108
-
109
-
echo " building oxcaml universe (5.2.0+ox switch)..."
110
-
jtw opam --switch=5.2.0+ox -o "$UNIVERSES/oxcaml"
111
-
112
-
for d in universe universe-v2 universe-v3 universe-oxcaml; do
113
-
rm -rf "$DEMO_DIR/$d"
114
-
done
115
-
116
-
cp -r "$UNIVERSES/default" "$DEMO_DIR/universe"
117
-
echo " deployed universe/"
118
-
119
-
cp -r "$UNIVERSES/default" "$DEMO_DIR/universe-v2"
120
-
echo " deployed universe-v2/"
121
-
122
-
cp -r "$UNIVERSES/v3" "$DEMO_DIR/universe-v3"
123
-
echo " deployed universe-v3/"
124
-
125
-
cp -r "$UNIVERSES/oxcaml" "$DEMO_DIR/universe-oxcaml"
126
-
echo " deployed universe-oxcaml/"
127
-
128
-
for d in universe universe-v2 universe-v3 universe-oxcaml; do
129
-
cp "$SITE/_x-ocaml/x-ocaml.js" "$DEMO_DIR/$d/x-ocaml.js"
130
-
done
131
-
fi
132
-
133
-
echo ""
134
-
echo "=== Done ==="
135
-
echo ""
136
-
echo "Site root: $SITE/"
137
-
echo ""
138
-
echo "Key pages:"
139
-
echo " /index.html — site home"
140
-
echo " /blog/index.html — blog index"
141
-
echo " /notebooks/index.html — notebooks index"
142
-
echo " /notebooks/foundations/index.html — foundations of CS"
143
-
echo " /projects/index.html — projects"
144
-
echo " /reference/ — API reference docs"
145
-
echo " /reference/odoc-interactive-extension/ — interactive demos"
146
-
147
-
if $SERVE; then
148
-
echo ""
149
-
echo "Starting HTTP server on http://localhost:8080"
150
-
cd "$SITE"
151
-
exec python3 -m http.server 8080
152
-
fi
30
+
exec "$MONO/build-site.sh" "${BUILD_ARGS[@]}"
+7
-12
js_top_worker/bin/jtw.ml
+7
-12
js_top_worker/bin/jtw.ml
···
360
360
List.iter
361
361
(fun (dir, dcs) ->
362
362
let findlib_dir = Ocamlfind.findlib_dir () |> Fpath.v in
363
-
let d = Fpath.relativize ~root:findlib_dir dir in
364
-
match d with
365
-
| None ->
366
-
Format.eprintf "Failed to relativize %a wrt %a\n%!" Fpath.pp dir
367
-
Fpath.pp findlib_dir
368
-
| Some dir ->
369
-
Format.eprintf "Generating %a\n%!" Fpath.pp dir;
370
-
let dir = Fpath.(output_dir / "lib" // dir) in
371
-
let _ = Bos.OS.Dir.create dir in
372
-
let oc = open_out Fpath.(dir / "dynamic_cmis.json" |> to_string) in
373
-
Printf.fprintf oc "%s" dcs;
374
-
close_out oc)
363
+
let d = relativize_or_fallback ~findlib_dir dir in
364
+
Format.eprintf "Generating %a\n%!" Fpath.pp d;
365
+
let dir = Fpath.(output_dir / "lib" // d) in
366
+
let _ = Bos.OS.Dir.create dir in
367
+
let oc = open_out Fpath.(dir / "dynamic_cmis.json" |> to_string) in
368
+
Printf.fprintf oc "%s" dcs;
369
+
close_out oc)
375
370
init_cmis;
376
371
Format.eprintf "Number of cmis: %d\n%!" (List.length init_cmis);
377
372
+15
-15
site/notebooks/interactive_map.mld
+15
-15
site/notebooks/interactive_map.mld
···
29
29
(* Shared state *)
30
30
let bbox : Geotessera.bbox option ref = ref None
31
31
let mosaic : (Linalg.mat * int * int) option ref = ref None
32
+
let mosaic_bounds : Geotessera.bbox option ref = ref None
32
33
let projected : Linalg.mat option ref = ref None
33
34
let training_points : (float * float * int) list ref = ref []
34
35
let current_class = ref 0
···
119
120
| None -> Widget.update ~id:"status" (status_view "Error: draw a bounding box first!")
120
121
| Some b ->
121
122
Widget.update ~id:"status" (status_view "Fetching embeddings...");
122
-
let mat_full, h_full, w_full = Geotessera_jsoo.fetch_mosaic ~year:2024 b in
123
+
let mat_full, h_full, w_full, mosaic_bbox = Geotessera_jsoo.fetch_mosaic ~year:2024 b in
123
124
Widget.update ~id:"status"
124
125
(status_view (Printf.sprintf "Fetched %d×%d mosaic. Downsampling..." h_full w_full));
125
126
let mat, h, w = downsample_mosaic mat_full ~h:h_full ~w:w_full ~max_pixels:50_000 in
126
127
mosaic := Some (mat, h, w);
127
-
(* Compute actual mosaic bounds from tile grid (not user bbox) *)
128
-
let snap = Geotessera.snap_to_grid in
129
-
let mosaic_south = snap b.min_lat -. 0.05 in
130
-
let mosaic_north = snap b.max_lat +. 0.05 in
131
-
let mosaic_west = snap b.min_lon -. 0.05 in
132
-
let mosaic_east = snap b.max_lon +. 0.05 in
128
+
mosaic_bounds := Some mosaic_bbox;
129
+
(* Use actual WGS84 bounds computed from UTM tile geometry *)
130
+
let mosaic_south = mosaic_bbox.Geotessera.min_lat in
131
+
let mosaic_north = mosaic_bbox.Geotessera.max_lat in
132
+
let mosaic_west = mosaic_bbox.Geotessera.min_lon in
133
+
let mosaic_east = mosaic_bbox.Geotessera.max_lon in
133
134
Widget.update ~id:"status"
134
135
(status_view (Printf.sprintf "Working at %d×%d (%d pixels). Computing PCA..." h w (h * w)));
135
136
(* PCA to 3 components for RGB visualisation *)
···
228
229
229
230
{@ocaml x[
230
231
let () =
231
-
match !mosaic, !projected, !bbox with
232
-
| Some (mat, h, w), Some _proj, Some b ->
232
+
match !mosaic, !projected, !mosaic_bounds with
233
+
| Some (mat, h, w), Some _proj, Some mb ->
233
234
let points = !training_points in
234
235
if points = [] then
235
236
Widget.update ~id:"status" (status_view "Error: place some training points first!")
236
237
else begin
237
238
Widget.update ~id:"status"
238
239
(status_view (Printf.sprintf "Classifying with %d training points..." (List.length points)));
239
-
(* Compute actual mosaic bounds from tile grid *)
240
-
let snap = Geotessera.snap_to_grid in
241
-
let mosaic_south = snap b.min_lat -. 0.05 in
242
-
let mosaic_north = snap b.max_lat +. 0.05 in
243
-
let mosaic_west = snap b.min_lon -. 0.05 in
244
-
let mosaic_east = snap b.max_lon +. 0.05 in
240
+
(* Use stored mosaic bounds computed from UTM tile geometry *)
241
+
let mosaic_south = mb.Geotessera.min_lat in
242
+
let mosaic_north = mb.Geotessera.max_lat in
243
+
let mosaic_west = mb.Geotessera.min_lon in
244
+
let mosaic_east = mb.Geotessera.max_lon in
245
245
(* Convert geo coords to pixel coords and extract embeddings *)
246
246
let n_train = List.length points in
247
247
let train_mat = Linalg.create_mat ~rows:n_train ~cols:mat.Linalg.cols in
+2
-1
tessera-geotessera-jsoo/lib/geotessera_jsoo.mli
+2
-1
tessera-geotessera-jsoo/lib/geotessera_jsoo.mli
···
9
9
@raise Failure on non-200 status or if response cannot be read. *)
10
10
11
11
val fetch_mosaic : ?base_url:string -> ?version:string -> year:int ->
12
-
Geotessera.bbox -> Linalg.mat * int * int
12
+
Geotessera.bbox -> Linalg.mat * int * int * Geotessera.bbox
13
13
(** Fetch and assemble tiles for a bounding box using browser XHR.
14
+
Returns [(mosaic_mat, height, width, mosaic_bounds)].
14
15
Convenience wrapper around {!Geotessera.fetch_mosaic_sync}
15
16
with {!fetch} as the fetch function. *)
+96
-13
tessera-geotessera/lib/geotessera.ml
+96
-13
tessera-geotessera/lib/geotessera.ml
···
82
82
done;
83
83
result
84
84
85
+
(* Tile bounds and reprojection *)
86
+
87
+
let tile_bounds_wgs84 coord ~tile_h:_ ~tile_w:_ =
88
+
(* Nominal WGS84 bounds of the tile. The actual pixel data is on a UTM grid
89
+
that extends slightly beyond this, but after reprojection the output
90
+
is resampled onto this exact WGS84 grid. *)
91
+
{ min_lon = coord.lon -. 0.05;
92
+
min_lat = coord.lat -. 0.05;
93
+
max_lon = coord.lon +. 0.05;
94
+
max_lat = coord.lat +. 0.05 }
95
+
96
+
(* Reproject a single tile from its native UTM grid onto a regular WGS84 grid.
97
+
The input tile has [tile_h * tile_w] pixels on a UTM grid. We resample onto
98
+
[out_h * out_w] pixels covering the nominal 0.1-degree WGS84 tile box.
99
+
Uses nearest-neighbor interpolation. *)
100
+
let reproject_tile coord ~tile_h ~tile_w ~n_features ~out_h ~out_w data =
101
+
let zone = Utm.zone_of_lon coord.lon in
102
+
(* Compute the UTM bounding box of the 0.1-degree WGS84 tile *)
103
+
let corners = [
104
+
(coord.lon -. 0.05, coord.lat -. 0.05);
105
+
(coord.lon +. 0.05, coord.lat -. 0.05);
106
+
(coord.lon -. 0.05, coord.lat +. 0.05);
107
+
(coord.lon +. 0.05, coord.lat +. 0.05);
108
+
] in
109
+
let utm_corners = List.map (fun (lon, lat) -> Utm.wgs84_to_utm ~zone lon lat) corners in
110
+
let utm_min_e = ref Float.max_float in
111
+
let utm_max_e = ref Float.neg_infinity in
112
+
let utm_min_n = ref Float.max_float in
113
+
let utm_max_n = ref Float.neg_infinity in
114
+
List.iter (fun (e, n) ->
115
+
if e < !utm_min_e then utm_min_e := e;
116
+
if e > !utm_max_e then utm_max_e := e;
117
+
if n < !utm_min_n then utm_min_n := n;
118
+
if n > !utm_max_n then utm_max_n := n;
119
+
) utm_corners;
120
+
(* UTM pixel spacing in the input tile *)
121
+
let utm_pixel_e = (!utm_max_e -. !utm_min_e) /. Float.of_int tile_w in
122
+
let utm_pixel_n = (!utm_max_n -. !utm_min_n) /. Float.of_int tile_h in
123
+
(* Output WGS84 grid *)
124
+
let wgs_south = coord.lat -. 0.05 in
125
+
let wgs_north = coord.lat +. 0.05 in
126
+
let wgs_west = coord.lon -. 0.05 in
127
+
let wgs_east = coord.lon +. 0.05 in
128
+
let result = Linalg.create_mat ~rows:(out_h * out_w) ~cols:n_features in
129
+
for oi = 0 to out_h - 1 do
130
+
(* Latitude of this output row (row 0 = north) *)
131
+
let lat = wgs_north -. (Float.of_int oi +. 0.5) *. (wgs_north -. wgs_south) /. Float.of_int out_h in
132
+
for oj = 0 to out_w - 1 do
133
+
(* Longitude of this output column *)
134
+
let lon = wgs_west +. (Float.of_int oj +. 0.5) *. (wgs_east -. wgs_west) /. Float.of_int out_w in
135
+
(* Convert to UTM *)
136
+
let (e, n) = Utm.wgs84_to_utm ~zone lon lat in
137
+
(* Find nearest input pixel.
138
+
GeoTessera tiles have row 0 = north (northing decreases with row),
139
+
matching the landmask GeoTIFF Affine transform convention. *)
140
+
let pi = Float.to_int (Float.floor ((!utm_max_n -. n) /. utm_pixel_n)) in
141
+
let pj = Float.to_int (Float.floor ((e -. !utm_min_e) /. utm_pixel_e)) in
142
+
if pi >= 0 && pi < tile_h && pj >= 0 && pj < tile_w then begin
143
+
let in_idx = pi * tile_w + pj in
144
+
let out_idx = oi * out_w + oj in
145
+
for f = 0 to n_features - 1 do
146
+
Linalg.mat_set result out_idx f (Linalg.mat_get data in_idx f)
147
+
done
148
+
end
149
+
done
150
+
done;
151
+
result
152
+
85
153
(* Mosaic *)
86
154
87
155
let mosaic tiles bbox =
88
156
match tiles with
89
-
| [] -> (Linalg.create_mat ~rows:0 ~cols:0, 0, 0)
157
+
| [] -> (Linalg.create_mat ~rows:0 ~cols:0, 0, 0, { min_lon = 0.0; min_lat = 0.0; max_lon = 0.0; max_lat = 0.0 })
90
158
| (_, first_data, tile_h, tile_w) :: _ ->
91
159
let n_features = first_data.Linalg.cols in
92
-
(* Find grid extent *)
160
+
let _ = bbox in
161
+
(* Find grid extent from tile coordinates *)
93
162
let min_lat = ref Float.max_float in
94
163
let max_lat = ref Float.neg_infinity in
95
164
let min_lon = ref Float.max_float in
96
165
let max_lon = ref Float.neg_infinity in
97
-
let _ = bbox in
98
166
List.iter (fun (coord, _, _, _) ->
99
167
if coord.lat < !min_lat then min_lat := coord.lat;
100
168
if coord.lat > !max_lat then max_lat := coord.lat;
···
103
171
) tiles;
104
172
let n_tile_rows = Float.to_int (Float.round ((!max_lat -. !min_lat) /. 0.1)) + 1 in
105
173
let n_tile_cols = Float.to_int (Float.round ((!max_lon -. !min_lon) /. 0.1)) + 1 in
106
-
let total_h = n_tile_rows * tile_h in
107
-
let total_w = n_tile_cols * tile_w in
174
+
(* Choose output grid size per tile: match the input resolution.
175
+
Use the first tile's dimensions as the output size for all tiles. *)
176
+
let out_h = tile_h in
177
+
let out_w = tile_w in
178
+
let total_h = n_tile_rows * out_h in
179
+
let total_w = n_tile_cols * out_w in
108
180
let output = Linalg.create_mat ~rows:(total_h * total_w) ~cols:n_features in
109
181
List.iter (fun (coord, data, th, tw) ->
110
-
(* lat increases north, but row 0 is northernmost *)
182
+
(* Reproject this tile from UTM to a regular WGS84 grid *)
183
+
let reprojected = reproject_tile coord ~tile_h:th ~tile_w:tw
184
+
~n_features ~out_h ~out_w data in
185
+
(* Place reprojected tile in the mosaic grid *)
111
186
let gi = Float.to_int (Float.round ((!max_lat -. coord.lat) /. 0.1)) in
112
187
let gj = Float.to_int (Float.round ((coord.lon -. !min_lon) /. 0.1)) in
113
-
for pi = 0 to th - 1 do
114
-
for pj = 0 to tw - 1 do
115
-
let out_row = (gi * tile_h + pi) * total_w + (gj * tile_w + pj) in
116
-
let in_row = pi * tw + pj in
188
+
for pi = 0 to out_h - 1 do
189
+
for pj = 0 to out_w - 1 do
190
+
let out_row = (gi * out_h + pi) * total_w + (gj * out_w + pj) in
191
+
let in_row = pi * out_w + pj in
117
192
for f = 0 to n_features - 1 do
118
-
Linalg.mat_set output out_row f (Linalg.mat_get data in_row f)
193
+
Linalg.mat_set output out_row f (Linalg.mat_get reprojected in_row f)
119
194
done
120
195
done
121
196
done
122
197
) tiles;
123
-
(output, total_h, total_w)
198
+
(* Mosaic bounds = nominal WGS84 grid *)
199
+
let mosaic_bbox = {
200
+
min_lon = !min_lon -. 0.05;
201
+
min_lat = !min_lat -. 0.05;
202
+
max_lon = !max_lon +. 0.05;
203
+
max_lat = !max_lat +. 0.05;
204
+
} in
205
+
(output, total_h, total_w, mosaic_bbox)
124
206
125
207
(* High-level API *)
126
208
···
143
225
let tile_w = shape.(1) in
144
226
(coord, mat, tile_h, tile_w)
145
227
) tiles in
146
-
mosaic tile_data bbox
228
+
let (mat, h, w, mosaic_bbox) = mosaic tile_data bbox in
229
+
(mat, h, w, mosaic_bbox)
+13
-3
tessera-geotessera/lib/geotessera.mli
+13
-3
tessera-geotessera/lib/geotessera.mli
···
36
36
(** [int8 * float32_scale -> float32]. Input embeddings shape (H,W,128), scales shape (H,W).
37
37
Output: mat with rows=H*W, cols=128. *)
38
38
39
+
(** {1 Tile bounds} *)
40
+
41
+
val tile_bounds_wgs84 : tile_coord -> tile_h:int -> tile_w:int -> bbox
42
+
(** Compute the actual WGS84 bounding box of a tile's pixel grid.
43
+
The pixel grid is the UTM bounding box of the WGS84 tile corners,
44
+
which is wider than the nominal 0.1-degree tile due to UTM projection. *)
45
+
39
46
(** {1 Mosaicing} *)
40
47
41
-
val mosaic : (tile_coord * Linalg.mat * int * int) list -> bbox -> Linalg.mat * int * int
48
+
val mosaic : (tile_coord * Linalg.mat * int * int) list -> bbox -> Linalg.mat * int * int * bbox
42
49
(** Assemble tiles into a single matrix. Each tile is [(coord, data, tile_h, tile_w)].
43
-
Returns [(mosaic_mat, total_h, total_w)]. *)
50
+
Returns [(mosaic_mat, total_h, total_w, mosaic_bbox)].
51
+
The fourth element is the WGS84 bounding box of the mosaic, computed as the
52
+
union of all tile bounds via {!tile_bounds_wgs84}. *)
44
53
45
54
(** {1 High-level API} *)
46
55
···
50
59
?version:string ->
51
60
year:int ->
52
61
bbox ->
53
-
Linalg.mat * int * int
62
+
Linalg.mat * int * int * bbox
54
63
(** Fetch all tiles for a bbox, dequantize, and mosaic.
55
64
[fetch] is a blocking function that retrieves URL contents as a string.
65
+
Returns [(mosaic_mat, total_h, total_w, mosaic_bbox)].
56
66
Default base_url: ["https://dl2.geotessera.org"].
57
67
Default version: ["v1"]. *)
+105
tessera-geotessera/lib/utm.ml
+105
tessera-geotessera/lib/utm.ml
···
1
+
(* UTM projection utilities — WGS84 ellipsoid, Northern hemisphere *)
2
+
3
+
let a = 6378137.0
4
+
let f = 1.0 /. 298.257223563
5
+
let b_ = a *. (1.0 -. f)
6
+
let e2 = 2.0 *. f -. f *. f
7
+
let e'2 = e2 /. (1.0 -. e2)
8
+
let n = (a -. b_) /. (a +. b_)
9
+
let k0 = 0.9996
10
+
11
+
let zone_of_lon lon =
12
+
Float.to_int (Float.floor ((lon +. 180.0) /. 6.0)) + 1
13
+
14
+
let central_meridian zone =
15
+
Float.of_int (zone * 6 - 183)
16
+
17
+
let wgs84_to_utm ~zone lon lat =
18
+
let lon0 = central_meridian zone in
19
+
let lat_r = lat *. Float.pi /. 180.0 in
20
+
let lon_r = lon *. Float.pi /. 180.0 in
21
+
let lon0_r = lon0 *. Float.pi /. 180.0 in
22
+
let cos_lat = Float.cos lat_r in
23
+
let sin_lat = Float.sin lat_r in
24
+
let tan_lat = Float.tan lat_r in
25
+
let n2 = n *. n in
26
+
let n3 = n2 *. n in
27
+
let n4 = n3 *. n in
28
+
let n5 = n4 *. n in
29
+
let _n6 = n5 *. n in
30
+
(* Meridional arc using series expansion *)
31
+
let a0 = 1.0 +. n2 /. 4.0 +. n4 /. 64.0 in
32
+
let a2 = 1.5 *. (n -. n3 /. 8.0 +. n5 /. 64.0) in
33
+
let a4 = 15.0 /. 16.0 *. (n2 -. n4 /. 4.0) in
34
+
let a6 = 35.0 /. 48.0 *. (n3 -. 5.0 *. n5 /. 16.0) in
35
+
let a8 = 315.0 /. 512.0 *. n4 in
36
+
let a10 = 693.0 /. 1280.0 *. n5 in
37
+
let _ = a10 in
38
+
let big_a = a /. (1.0 +. n) *. a0 in
39
+
let s = big_a *. (lat_r
40
+
-. a2 *. Float.sin (2.0 *. lat_r)
41
+
+. a4 *. Float.sin (4.0 *. lat_r)
42
+
-. a6 *. Float.sin (6.0 *. lat_r)
43
+
+. a8 *. Float.sin (8.0 *. lat_r)) in
44
+
let nu = a /. Float.sqrt (1.0 -. e2 *. sin_lat *. sin_lat) in
45
+
let t = tan_lat in
46
+
let t2 = t *. t in
47
+
let t4 = t2 *. t2 in
48
+
let c = e'2 *. cos_lat *. cos_lat in
49
+
let dl = lon_r -. lon0_r in
50
+
let dl2 = dl *. dl in
51
+
let dl3 = dl2 *. dl in
52
+
let dl4 = dl3 *. dl in
53
+
let dl5 = dl4 *. dl in
54
+
let dl6 = dl5 *. dl in
55
+
let easting = 500000.0 +. k0 *. nu *. cos_lat *. (dl
56
+
+. (1.0 -. t2 +. c) *. cos_lat *. cos_lat *. dl3 /. 6.0
57
+
+. (5.0 -. 18.0 *. t2 +. t4 +. 72.0 *. c -. 58.0 *. e'2)
58
+
*. cos_lat *. cos_lat *. cos_lat *. cos_lat *. dl5 /. 120.0) in
59
+
let northing = k0 *. (s +. nu *. tan_lat *. (
60
+
cos_lat *. cos_lat *. dl2 /. 2.0
61
+
+. (5.0 -. t2 +. 9.0 *. c +. 4.0 *. c *. c) *. cos_lat *. cos_lat *. cos_lat *. cos_lat *. dl4 /. 24.0
62
+
+. (61.0 -. 58.0 *. t2 +. t4 +. 600.0 *. c -. 330.0 *. e'2)
63
+
*. cos_lat *. cos_lat *. cos_lat *. cos_lat *. cos_lat *. cos_lat *. dl6 /. 720.0)) in
64
+
(easting, northing)
65
+
66
+
let utm_to_wgs84 ~zone easting northing =
67
+
let lon0 = central_meridian zone in
68
+
let lon0_r = lon0 *. Float.pi /. 180.0 in
69
+
let e1 = (1.0 -. Float.sqrt (1.0 -. e2)) /. (1.0 +. Float.sqrt (1.0 -. e2)) in
70
+
let e1_2 = e1 *. e1 in
71
+
let e1_3 = e1_2 *. e1 in
72
+
let e1_4 = e1_3 *. e1 in
73
+
let x = easting -. 500000.0 in
74
+
let m = northing /. k0 in
75
+
let mu = m /. (a *. (1.0 -. e2 /. 4.0 -. 3.0 *. e2 *. e2 /. 64.0 -. 5.0 *. e2 *. e2 *. e2 /. 256.0)) in
76
+
let phi1 = mu
77
+
+. (3.0 *. e1 /. 2.0 -. 27.0 *. e1_3 /. 32.0) *. Float.sin (2.0 *. mu)
78
+
+. (21.0 *. e1_2 /. 16.0 -. 55.0 *. e1_4 /. 32.0) *. Float.sin (4.0 *. mu)
79
+
+. (151.0 *. e1_3 /. 96.0) *. Float.sin (6.0 *. mu)
80
+
+. (1097.0 *. e1_4 /. 512.0) *. Float.sin (8.0 *. mu) in
81
+
let sin_phi1 = Float.sin phi1 in
82
+
let cos_phi1 = Float.cos phi1 in
83
+
let tan_phi1 = Float.tan phi1 in
84
+
let nu1 = a /. Float.sqrt (1.0 -. e2 *. sin_phi1 *. sin_phi1) in
85
+
let rho1 = a *. (1.0 -. e2) /. (Float.pow (1.0 -. e2 *. sin_phi1 *. sin_phi1) 1.5) in
86
+
let t1 = tan_phi1 in
87
+
let t1_2 = t1 *. t1 in
88
+
let c1 = e'2 *. cos_phi1 *. cos_phi1 in
89
+
let d = x /. (nu1 *. k0) in
90
+
let d2 = d *. d in
91
+
let d3 = d2 *. d in
92
+
let d4 = d3 *. d in
93
+
let d5 = d4 *. d in
94
+
let d6 = d5 *. d in
95
+
let lat = phi1 -. (nu1 *. tan_phi1 /. rho1) *. (
96
+
d2 /. 2.0
97
+
-. (5.0 +. 3.0 *. t1_2 +. 10.0 *. c1 -. 4.0 *. c1 *. c1 -. 9.0 *. e'2) *. d4 /. 24.0
98
+
+. (61.0 +. 90.0 *. t1_2 +. 298.0 *. c1 +. 45.0 *. t1_2 *. t1_2 -. 252.0 *. e'2 -. 3.0 *. c1 *. c1) *. d6 /. 720.0) in
99
+
let lon = lon0_r +. (d
100
+
-. (1.0 +. 2.0 *. t1_2 +. c1) *. d3 /. 6.0
101
+
+. (5.0 -. 2.0 *. c1 +. 28.0 *. t1_2 -. 3.0 *. c1 *. c1 +. 8.0 *. e'2 +. 24.0 *. t1_2 *. t1_2) *. d5 /. 120.0)
102
+
/. cos_phi1 in
103
+
let lat_deg = lat *. 180.0 /. Float.pi in
104
+
let lon_deg = lon *. 180.0 /. Float.pi in
105
+
(lon_deg, lat_deg)
+15
tessera-geotessera/lib/utm.mli
+15
tessera-geotessera/lib/utm.mli
···
1
+
(** UTM projection utilities — WGS84 ellipsoid, Northern hemisphere. *)
2
+
3
+
val zone_of_lon : float -> int
4
+
(** UTM zone number from longitude. *)
5
+
6
+
val central_meridian : int -> float
7
+
(** Central meridian (degrees) of a UTM zone. *)
8
+
9
+
val wgs84_to_utm : zone:int -> float -> float -> float * float
10
+
(** [wgs84_to_utm ~zone lon lat] converts WGS84 (lon, lat) in degrees
11
+
to UTM (easting, northing) in metres. Northern hemisphere only. *)
12
+
13
+
val utm_to_wgs84 : zone:int -> float -> float -> float * float
14
+
(** [utm_to_wgs84 ~zone easting northing] converts UTM (easting, northing)
15
+
in metres to WGS84 (lon, lat) in degrees. Northern hemisphere only. *)
+1
-4
tessera-geotessera/test/dune
+1
-4
tessera-geotessera/test/dune
+3
-3
tessera-geotessera/test/test_geotessera.ml
+3
-3
tessera-geotessera/test/test_geotessera.ml
···
119
119
(coord_a, tile_a, tile_h, tile_w);
120
120
(coord_b, tile_b, tile_h, tile_w);
121
121
] in
122
-
let (result, rh, rw) = Geotessera.mosaic tiles bbox in
122
+
let (result, rh, rw, _mosaic_bbox) = Geotessera.mosaic tiles bbox in
123
123
(* total_h = 1 tile row * 2 = 2, total_w = 2 tile cols * 3 = 6 *)
124
124
Alcotest.(check int) "total_h" 2 rh;
125
125
Alcotest.(check int) "total_w" 6 rw;
···
155
155
(coord_north, north_tile, tile_h, tile_w);
156
156
(coord_south, south_tile, tile_h, tile_w);
157
157
] in
158
-
let (result, rh, rw) = Geotessera.mosaic tiles bbox in
158
+
let (result, rh, rw, _mosaic_bbox) = Geotessera.mosaic tiles bbox in
159
159
Alcotest.(check int) "total_h" 4 rh;
160
160
Alcotest.(check int) "total_w" 2 rw;
161
161
(* North tile (higher lat) should be at row 0 (top) *)
···
181
181
min_lon = 0.12; min_lat = 52.12;
182
182
max_lon = 0.18; max_lat = 52.18;
183
183
} in
184
-
let (result, rh, rw) = Geotessera.fetch_mosaic_sync
184
+
let (result, rh, rw, _mosaic_bbox) = Geotessera.fetch_mosaic_sync
185
185
~fetch ~base_url:"http://mock" ~version:"v1" ~year:2021 bbox in
186
186
(* tile is 2x2 with 2 features *)
187
187
Alcotest.(check int) "h" 2 rh;
+1
-1
test-pipeline/test_pipeline.ml
+1
-1
test-pipeline/test_pipeline.ml
···
189
189
in
190
190
191
191
let bbox = Geotessera.{ min_lon = 0.0; min_lat = 52.0; max_lon = 0.05; max_lat = 52.05 } in
192
-
let mosaic, h, w = Geotessera.fetch_mosaic_sync ~fetch ~year:2024 bbox in
192
+
let mosaic, h, w, _mosaic_bbox = Geotessera.fetch_mosaic_sync ~fetch ~year:2024 bbox in
193
193
194
194
Printf.printf "Mosaic from fetch: %d rows x %d cols, %dx%d pixels\n" mosaic.rows mosaic.cols h w;
195
195
Alcotest.(check int) "mosaic cols" 128 mosaic.cols;