forked from
anil.recoil.org/monopam-myspace
My aggregated monorepo of OCaml code, automaintained
Fix Zarr reprojection to use actual region bounds
reproject_tile assumed a hardcoded 0.1° WGS84 window centred on the
bbox centre, regardless of the actual fetched region size. This
caused the overlay to be misaligned with the map.
Replace with inline reprojection that uses the actual UTM pixel
bounds and the derived WGS84 bounds. Each output pixel's WGS84
coordinate is converted to UTM and mapped to the nearest input
pixel using the real pixel grid, not a fixed tile assumption.
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
+34
-9
+34
-9
tessera-zarr/lib/tessera_zarr.ml
+34
-9
tessera-zarr/lib/tessera_zarr.ml
···
36
36
let open Geotessera in
37
37
(* 1. Determine UTM zone from bbox centre *)
38
38
let center_lon = (bbox.min_lon +. bbox.max_lon) /. 2.0 in
39
-
let center_lat = (bbox.min_lat +. bbox.max_lat) /. 2.0 in
40
39
let zone = Geotessera.Utm.zone_of_lon center_lon in
41
40
let zone_name = Printf.sprintf "utm%d" zone in
42
41
···
109
108
done
110
109
done;
111
110
112
-
(* 6. Reproject from UTM to WGS84 *)
113
-
let out_h = tile_h in
114
-
let out_w = tile_w in
115
-
let coord = { lon = center_lon; lat = center_lat } in
116
-
let reprojected = Geotessera.reproject_tile coord
117
-
~tile_h ~tile_w ~n_features ~out_h ~out_w mat in
118
-
119
-
(* 7. Compute WGS84 bounds *)
111
+
(* 6. Compute actual UTM bounds of fetched pixels *)
120
112
let utm_w = zi.origin_easting +. Float.of_int col_start *. pixel_size in
121
113
let utm_e = zi.origin_easting +. Float.of_int col_stop *. pixel_size in
122
114
let utm_n = zi.origin_northing -. Float.of_int row_start *. pixel_size in
123
115
let utm_s = zi.origin_northing -. Float.of_int row_stop *. pixel_size in
116
+
117
+
(* 7. Compute WGS84 bounds from the UTM corners *)
124
118
let (w_lon, s_lat) = Geotessera.Utm.utm_to_wgs84 ~zone utm_w utm_s in
125
119
let (e_lon, n_lat) = Geotessera.Utm.utm_to_wgs84 ~zone utm_e utm_n in
126
120
let wgs_bbox = { min_lon = w_lon; min_lat = s_lat;
127
121
max_lon = e_lon; max_lat = n_lat } in
122
+
123
+
(* 8. Reproject from UTM to WGS84 using the actual bounds.
124
+
For each pixel in the output WGS84 grid, convert to UTM and
125
+
sample the nearest input pixel. *)
126
+
progress "Reprojecting...";
127
+
let out_h = tile_h in
128
+
let out_w = tile_w in
129
+
let reprojected = Linalg.create_mat ~rows:(out_h * out_w) ~cols:n_features in
130
+
for oi = 0 to out_h - 1 do
131
+
let lat = wgs_bbox.max_lat -.
132
+
(Float.of_int oi +. 0.5) *. (wgs_bbox.max_lat -. wgs_bbox.min_lat)
133
+
/. Float.of_int out_h in
134
+
for oj = 0 to out_w - 1 do
135
+
let lon = wgs_bbox.min_lon +.
136
+
(Float.of_int oj +. 0.5) *. (wgs_bbox.max_lon -. wgs_bbox.min_lon)
137
+
/. Float.of_int out_w in
138
+
let (e, n) = Geotessera.Utm.wgs84_to_utm ~zone lon lat in
139
+
(* Map UTM coord to input pixel *)
140
+
let pi = Float.to_int (Float.floor ((utm_n -. n) /. pixel_size))
141
+
- row_start in
142
+
let pj = Float.to_int (Float.floor ((e -. utm_w) /. pixel_size)) in
143
+
if pi >= 0 && pi < tile_h && pj >= 0 && pj < tile_w then begin
144
+
let in_idx = pi * tile_w + pj in
145
+
let out_idx = oi * out_w + oj in
146
+
for f = 0 to n_features - 1 do
147
+
Linalg.mat_set reprojected out_idx f
148
+
(Linalg.mat_get mat in_idx f)
149
+
done
150
+
end
151
+
done
152
+
done;
128
153
129
154
Lwt.return (reprojected, out_h, out_w, wgs_bbox)