Geotessera library for OCaml

Return transformation matrix too

Changed files
+29 -11
example
src
+3 -3
example/main.ml
···
let belfast =
Gt.Bbox.v ~min_lat:54.48 ~min_lon:(-6.11) ~max_lat:54.66 ~max_lon:(-5.78)
-
let pp_result fmt (crs, (point : Gt.point), emb) =
-
Fmt.pf fmt "Embedding at (%.2f, %.2f) with CRS:%i: %a" point.lat point.lon crs
-
Nx.pp_shape (Nx.shape emb)
+
let pp_result fmt ((crs, transform), (point : Gt.point), emb) =
+
Fmt.pf fmt "Embedding at (%.2f, %.2f) with CRS:%i: %a\ntransform: [%a]"
+
point.lat point.lon crs Nx.pp_shape (Nx.shape emb) Nx.pp transform
let () =
Mirage_crypto_rng_unix.use_default ();
+17 -5
src/geotessera.ml
···
([], [], []) paths
|> fun (ps, s, t) -> List.combine ps (List.combine s t)
-
let crs_of_landmask lm =
+
let transformation_matrix e =
+
let scale = Tiff.Ifd.pixel_scale e in
+
let tiepoint = Tiff.Ifd.tiepoint e in
+
let arr = Nx.zeros Nx.float64 [| 6 |] in
+
(* X *)
+
Nx.set_item [ 0 ] tiepoint.(3) arr;
+
(* sX *) Nx.set_item [ 1 ] scale.(0) arr;
+
Nx.set_item [ 3 ] tiepoint.(4) arr;
+
Nx.set_item [ 5 ] (-.scale.(1)) arr;
+
arr
+
+
let crs_and_transform_of_landmask lm =
Eio.Path.with_open_in lm @@ fun r ->
let ro = Eio.File.pread_exact r in
let tiff = Tiff.from_file Tiff.Float32 ro in
-
let geos = Tiff.ifd tiff |> Tiff.Ifd.geo_key_directory in
-
Tiff.Ifd.GeoKeys.projected_crs geos
+
let ifd = Tiff.ifd tiff in
+
let geos = Tiff.Ifd.geo_key_directory ifd in
+
(Tiff.Ifd.GeoKeys.projected_crs geos, transformation_matrix ifd)
let download_landmasks t v =
Eio.Switch.run @@ fun sw ->
···
(* Eio.traceln "%a\n%a" Digestif.SHA256.pp disk_hash *)
(* Digestif.SHA256.pp hash; *)
(* assert (Digestif.SHA256.equal disk_hash hash); *)
-
let crs = crs_of_landmask path in
+
let crs = crs_and_transform_of_landmask path in
(crs, path)
| `Delete_and_download | `Download ->
let () =
···
Client.with_body ~default:() t.client display ~uri
@@ fun body -> Flow.copy body w
in
-
let crs = crs_of_landmask path in
+
let crs = crs_and_transform_of_landmask path in
(crs, path)
| `Check_hash _ -> Fmt.failwith "%s exists but is not a file!" name)
v
+9 -3
src/geotessera.mli
···
_ env ->
year:int ->
Bbox.t ->
-
(int * point * Eio.Fs.dir_ty Eio.Path.t * (embedding * scales)) list
+
((int * (float, Bigarray.float64_elt) Nx.t)
+
* point
+
* Eio.Fs.dir_ty Eio.Path.t
+
* (embedding * scales))
+
list
(** [fetch env ~year bbox] returns a list of
-
[(crs, centre, landmask_path, scaled_embeddings)] for the tiles available in
-
[bbox] for [year]. *)
+
[((crs, transform), centre, landmask_path, scaled_embeddings)] for the tiles
+
available in [bbox] for [year].
+
+
Note that the [transformation] is returned in GDAL format. *)
val scale : embedding * scales -> (float, Bigarray.float32_elt) Nx.t
(** Scale up the embeddings *)