Cloud optimised geotiff main

CHANGES.md

@@ -27,6 +27,8 @@
 - Fix app crash when rendering feature info with a custom title.
 - Added new `CkanCatalogGroup` traits `resourceIdTemplate` and `restrictResourceIdTemplateToOrgsWithNames` to generate custom resource IDs for CKAN resources with unstable IDs.
 - Fix `acessType` resolution for `MagdaReference` so that it uses the default terria resolution strategy when `magdaRecord` is not defined.
+- Implement basic Cloud Optimised Geotiff support
+- [The next improvement]
 
 #### 8.2.28 - 2023-04-28
 

architecture/0011-leaflet-layer-override.md

@@ -0,0 +1,22 @@
+# 11. Leaflet Layer Override
+
+Date: 2023-06-26
+
+## Status
+
+Proposed
+
+## Context
+
+In order to support Cloud Optimised Geotiffs, we have encountered a situation where a Cesium Imagery Provider cannot provide imagery for both the Cesium viewer and the Leaflet viewer. We succesfully implemented display of COGs on the Cesium canvas with a library [TIFFImageryProvider](https://github.com/hongfaqiu/TIFFImageryProvider). However, this cannot be efficiently adapted to provide tiles for the Leaflet viewer. This is because Cesium is flexible with tiling schemes, but Leaflet is not.
+
+## Decision
+
+All mappable items now have an optional property to override the Imagery Provider when creating a Leaflet layer. This is currently `ImageryParts.overrideCreateLeafletLayer` and should return a `TerriaLeafletLayer` - a Leaflet layer with the required Terria enhancements such as Splitter support, Feature Picking support etc.
+
+When this property is specified, the Imagery Provider for the map item will still be set up when adding an item, and used in 3D mode to draw to the canvas. But the new ovveride function can specify a different way to draw to the Leaflet canvas. In the case of Cloud Optimised Geotiff catalog items, the overide function says that we should use the [Georaster Layer for Leaflet](https://github.com/GeoTIFF/georaster-layer-for-leaflet) library to create a Leaflet layer.
+
+## Consequences
+
+- This is less efficient than using the Imaery Provider to make requests for both 2D and 3D mode. In some situations, data will be re-requested if a user is switching between modes.
+- It is also more difficult for code maintenance, as functionality like Feature Picking must be implemented twice.

lib/Core/colourMappings.ts

@@ -0,0 +1,113 @@
+/** A collection of functions to map an array of band values at a pixel to an rgb colour */
+
+import chroma from "chroma-js";
+
+// Simply return Hex colour string from RGB. A useful starting point to modify.
+export function RGBAToHex(r: number, g: number, b: number) {
+  let rString = r.toString(16);
+  let gString = g.toString(16);
+  let bString = b.toString(16);
+
+  if (rString.length == 1) rString = "0" + rString;
+  if (gString.length == 1) gString = "0" + gString;
+  if (bString.length == 1) bString = "0" + bString;
+
+  return "#" + rString + gString + bString;
+}
+
+export function RBGToYCbCr(values: number[]) {
+  const r = Math.round(values[0] + 1.402 * (values[2] - 0x80));
+  const g = Math.round(
+    values[0] - 0.34414 * (values[1] - 0x80) - 0.71414 * (values[2] - 0x80)
+  );
+  const b = Math.round(values[0] + 1.772 * (values[1] - 0x80));
+  return `rgb(${r},${g},${b})`;
+}
+
+export function mapElevationToRgbaSmoothed(
+  values: number[],
+  waterLevel: number
+) {
+  // Graduated blues and oranges
+  // More examples of colour functions here: https://github.com/GeoTIFF/georaster-layer-for-leaflet/issues/106
+  const elevation = values[0];
+
+  const min = -5; // Hardcode min and max, TODO: get from the current part of the returned COG!
+  const max = 5;
+  const range = 10;
+
+  if (elevation < -50) {
+    return "";
+  } else if (elevation >= waterLevel) {
+    const scale = chroma.scale("oranges").domain([0, 1]);
+    const scaledPixelValue = (max - elevation) / range;
+    const color = scale(scaledPixelValue).hex();
+    return color;
+  } else if (elevation < waterLevel) {
+    const scale = chroma.scale("blues").domain([1, 0]);
+    const scaledPixelValue = (elevation - min) / range;
+    const color = scale(scaledPixelValue).hex();
+    return color;
+  } else {
+    return "";
+  }
+}
+
+export function rgbFromSeparateBands(red: number, green: number, blue: number) {
+  // Trying to map values from Int16 to 0-255...
+  // const max = Math.pow(2, 14);
+  const max = Math.pow(2, 14) / 16; // TODO: Had to divie by 16 to make it look right...  Why?
+  red = Math.round((255 * red) / max);
+  green = Math.round((255 * green) / max);
+  blue = Math.round((255 * blue) / max);
+
+  // make sure no values exceed 255
+  red = Math.min(red, 255);
+  green = Math.min(green, 255);
+  blue = Math.min(blue, 255);
+
+  // treat all black as no data
+  if (red === 0 && green === 0 && blue === 0) return "";
+  return `rgb(${red}, ${green}, ${blue})`;
+}
+
+export function trueColour(red: number, green: number, blue: number) {
+  // Trying to map values from Int16 to 0-255...
+  // const max = Math.pow(2, 14);
+  const max = Math.pow(2, 14) / 16; // TODO: Had to divie by 16 to make it look right...  Why?
+  red = Math.round((255 * red) / max);
+  green = Math.round((255 * green) / max);
+  blue = Math.round((255 * blue) / max);
+
+  // make sure no values exceed 255
+  red = Math.min(red, 255);
+  green = Math.min(green, 255);
+  blue = Math.min(blue, 255);
+
+  // treat all black as no data
+  if (red === 0 && green === 0 && blue === 0) return "";
+  return `rgb(${red}, ${green}, ${blue})`;
+}
+
+const scale = chroma
+  .scale(["#C7BB95", "#FEFEE1", "#6E9F62", "#032816", "black"])
+  .domain([0, 0.2, 0.4, 0.6, 0.8]);
+
+export function ndvi(red: number, nir: number) {
+  const dividend = nir - red;
+  const divisor = nir + red;
+  let result;
+  if (dividend === 0 && divisor === 0) {
+    // probably no reading here
+    return "rgba(0,0,0,0)";
+  }
+  if (dividend === 0 || divisor === 0) {
+    result = 0;
+  } else {
+    result = dividend === 0 ? 0 : dividend / divisor;
+  }
+  if (result <= 0.1) return "blue";
+  if (result >= 0.8) return "black";
+
+  return scale(result).hex();
+}

lib/Core/getDataType.ts

@@ -124,6 +124,14 @@ const builtinRemoteDataTypes: RemoteDataType[] = [
   {
     value: "json",
     name: "core.dataType.json"
+  },
+  {
+    value: "cog",
+    name: "core.dataType.cog"
+  },
+  {
+    value: "cog-composite",
+    name: "core.dataType.cog-composite"
   }
   // Add next builtin remote upload type
 ];

lib/Map/Leaflet/GeorasterTerriaLayer.ts

@@ -0,0 +1,235 @@
+import L from "leaflet";
+import {
+  autorun,
+  computed,
+  IReactionDisposer,
+  makeObservable,
+  observable
+} from "mobx";
+import CesiumEvent from "terriajs-cesium/Source/Core/Event";
+import SplitDirection from "terriajs-cesium/Source/Scene/SplitDirection";
+import GeoRasterLayer, {
+  GeoRasterLayerOptions
+} from "georaster-layer-for-leaflet";
+import ImageryProvider from "terriajs-cesium/Source/Scene/ImageryProvider";
+import { identify } from "geoblaze";
+import ImageryLayerFeatureInfo from "terriajs-cesium/Source/Scene/ImageryLayerFeatureInfo";
+import Proj4Definitions from "../Vector/Proj4Definitions";
+import Leaflet from "../../Models/Leaflet";
+const proj4 = require("proj4").default;
+
+// TODO: Cannot extend GeoRasterLayerOptions why?
+// interface GeoRasterTerriaLayerOptions extends GeoRasterLayerOptions {
+//   imageryProvider: ImageryProvider;
+// }
+
+// We have to ts-ignore here as the type of GeoRasterLayer is expressed in an incompatible way
+//@ts-ignore
+
+export default class GeorasterTerriaLayer extends GeoRasterLayer {
+  readonly errorEvent: CesiumEvent = new CesiumEvent();
+  readonly initialized: boolean = false;
+  readonly _usable: boolean = false;
+  readonly _delayedUpdate: unknown = undefined;
+  readonly _zSubtract: number = 0;
+  readonly _previousCredits: unknown[] = [];
+
+  @observable splitDirection = SplitDirection.NONE;
+  @observable splitPosition: number = 0.5;
+
+  constructor(
+    private leaflet: Leaflet | undefined,
+    options: GeoRasterLayerOptions,
+    imageryProvider: ImageryProvider | undefined
+  ) {
+    super(Object.assign(options, { async: true, tileSize: 256 }));
+    this.imageryProvider = imageryProvider; // TODO: add to Options instead?
+
+    makeObservable(this);
+
+    // Handle splitter rection (and disposing reaction)
+    let disposeSplitterReaction: IReactionDisposer | undefined;
+    this.on("add", () => {
+      if (!disposeSplitterReaction) {
+        disposeSplitterReaction = this._reactToSplitterChange();
+      }
+    });
+    this.on("remove", () => {
+      if (disposeSplitterReaction) {
+        disposeSplitterReaction();
+        disposeSplitterReaction = undefined;
+      }
+    });
+  }
+
+  _reactToSplitterChange() {
+    return autorun(() => {
+      const container = this.getContainer();
+      if (container === null) {
+        return;
+      }
+
+      if (this.splitDirection === SplitDirection.LEFT) {
+        const { left: clipLeft } = this._clipsForSplitter;
+        container.style.clip = clipLeft;
+      } else if (this.splitDirection === SplitDirection.RIGHT) {
+        const { right: clipRight } = this._clipsForSplitter;
+        container.style.clip = clipRight;
+      } else {
+        container.style.clip = "auto";
+      }
+    });
+  }
+
+  @computed
+  get _clipsForSplitter() {
+    let clipLeft = "";
+    let clipRight = "";
+    let clipPositionWithinMap;
+    let clipX;
+
+    // TODO: Error in here - some undefined stuff
+    if (this.leaflet?.size && this.leaflet.nw && this.leaflet.se) {
+      clipPositionWithinMap = this.leaflet.size.x * this.splitPosition;
+      clipX = Math.round(this.leaflet.nw.x + clipPositionWithinMap);
+      clipLeft =
+        "rect(" +
+        [this.leaflet.nw.y, clipX, this.leaflet.se.y, this.leaflet.nw.x].join(
+          "px,"
+        ) +
+        "px)";
+      clipRight =
+        "rect(" +
+        [this.leaflet.nw.y, this.leaflet.se.x, this.leaflet.se.y, clipX].join(
+          "px,"
+        ) +
+        "px)";
+    }
+
+    return {
+      left: clipLeft,
+      right: clipRight,
+      clipPositionWithinMap: clipPositionWithinMap,
+      clipX: clipX
+    };
+  }
+
+  // Transform the feature picking coordinates to the native projection of the Georaster layer.
+  getFeaturePickingCoords(
+    map: L.Map,
+    longitudeRadians: number,
+    latitudeRadians: number
+  ) {
+    // get Georaster projection
+    const projection = this.extent.srs;
+    const reprojectFn = this.reprojectToSourceFn(projection);
+
+    // convert long and lat radians to x and y in native projection
+    const lat = (latitudeRadians / Math.PI) * 180;
+    const lon = (longitudeRadians / Math.PI) * 180;
+    const nativeCoords = reprojectFn([lon, lat]);
+
+    // Also include zoom level to fit previous implmentations
+    const level = Math.round(map.getZoom());
+
+    return {
+      x: nativeCoords[0],
+      y: nativeCoords[1],
+      level: level
+    };
+  }
+
+  reprojectToSourceFn = (sourceEpsgCode: string) => {
+    const sourceDef =
+      sourceEpsgCode in Proj4Definitions
+        ? new proj4.Proj(Proj4Definitions[sourceEpsgCode])
+        : undefined;
+
+    return proj4("EPSG:4326", sourceDef).forward;
+  };
+
+  /** Use the Geoblaze library to get pixel values by operating on the GeoRaster object.
+   * TODO: Is this giving the display values at that point, or the raw values of the full resolution raster at those coordinates?
+   * This is discussed in https://github.com/GeoTIFF/georaster-layer-for-leaflet/issues/104
+   * Currently his function is costly - `geoblaze.identify` takes time and appears to download the highest resolution tile for the area clicked.
+   * This is probably the only way if we want to get the actual raw pixel values at that point.
+   **/
+
+  async pickFeatures(
+    x: number,
+    y: number,
+    level: number,
+    longitudeRadians: number,
+    latitudeRadians: number
+  ) {
+    const featureInfo = new ImageryLayerFeatureInfo();
+    featureInfo.name = `lon:${((longitudeRadians / Math.PI) * 180).toFixed(
+      6
+    )}, lat:${((latitudeRadians / Math.PI) * 180).toFixed(6)}`;
+    const data: { [index: number]: any } = {};
+
+    for (let i = 0; i < this.georasters.length; i++) {
+      const res = await identify(this.georasters[i], [x, y]);
+
+      if (res) {
+        res.forEach((item: any): void => {
+          data[i] = item;
+        });
+        featureInfo.configureDescriptionFromProperties(data);
+      }
+    }
+
+    featureInfo.data = data;
+    return [featureInfo];
+  }
+
+  /** SUPERSEDED: These functions are alternative that do the feature picking using the TIFFImageryProvider.
+   * This is what the 3D Cesium map uses, but in the case of COG layers we are using a different layer for Leaflet 2D mode.
+   * It is more efficient to directly query the raster data fetched by Georaster Layer for Leaflet, that we are using on the 2D side.
+   * It can also work to query the pixel values using TIFFImageryProvider, but this requires downloading the tiles using that provider, and is inefficient.
+   * These functions have been retained during the testing phase of this new functionality, in case they offer other alternative.
+   */
+
+  //   // will get the coords in the tiling scheme provided by the TIFFImageryProvider
+  //   getFeaturePickingCoords(
+  //     map: L.Map,
+  //     longitudeRadians: number,
+  //     latitudeRadians: number
+  //   ) {
+  //     const ll = new Cartographic(
+  //       CesiumMath.negativePiToPi(longitudeRadians),
+  //       latitudeRadians,
+  //       0.0
+  //     );
+  //     const level = Math.round(map.getZoom());
+
+  //     return this.imageryProvider.readyPromise.then(() => {
+  //       const tilingScheme = this.imageryProvider.tilingScheme;
+  //       const coords = tilingScheme.positionToTileXY(ll, level);
+  //       return {
+  //         x: coords.x,
+  //         y: coords.y,
+  //         level: level
+  //       };
+  //     });
+  //   }
+
+  //   // THIS VERSION USES THE IMAGERY PROVIDER to get the values
+  //   pickFeatures(
+  //     x: number,
+  //     y: number,
+  //     level: number,
+  //     longitudeRadians: number,
+  //     latitudeRadians: number
+  //   ) {
+  //     return this.imageryProvider.readyPromise.then(() => {
+  //       return this.imageryProvider.pickFeatures(
+  //         x,
+  //         y,
+  //         level,
+  //         longitudeRadians,
+  //         latitudeRadians
+  //       );
+  //     });
+  //   }
+}