Reprojection with proj4 (#5)

* wip: implement reprojection on proj4

* update lockfile

* Finish off reprojection

* Add reprojection test

package.json

@@ -38,6 +38,7 @@
     "@rollup/plugin-terser": "^0.4.3",
     "@rollup/plugin-typescript": "^11.1.2",
     "@types/node": "^20.9.3",
+    "@types/proj4": "^2",
     "apache-arrow": "^14",
     "prettier": "^3.1.0",
     "rollup": "^4.1.5",
@@ -51,6 +52,7 @@
     "yarn": "4.0.2"
   },
   "dependencies": {
-    "@math.gl/polygon": "^4.0.0"
+    "@math.gl/polygon": "^4.0.0",
+    "proj4": "^2.9.2"
   }
 }

src/algorithm/coords.ts

@@ -0,0 +1,183 @@
+import * as arrow from "apache-arrow";
+import {
+  LineString,
+  MultiLineString,
+  MultiPoint,
+  MultiPolygon,
+  Point,
+  Polygon,
+} from "../type";
+import {
+  GeoArrowData,
+  LineStringData,
+  MultiLineStringData,
+  MultiPointData,
+  MultiPolygonData,
+  PointData,
+  PolygonData,
+  isLineStringData,
+  isMultiLineStringData,
+  isMultiPointData,
+  isMultiPolygonData,
+  isPointData,
+  isPolygonData,
+} from "../data";
+import {
+  getLineStringChild,
+  getMultiPolygonChild,
+  getPointChild,
+  getPolygonChild,
+} from "../child";
+import { assert, assertFalse } from "./utils/assert";
+
+// For now, simplify our lives by focusing on 2D
+type MapCoordsCallback = (x: number, y: number) => [number, number];
+
+export function mapCoords(
+  input: PointData,
+  callback: MapCoordsCallback,
+): PointData;
+export function mapCoords(
+  input: LineStringData,
+  callback: MapCoordsCallback,
+): LineStringData;
+export function mapCoords(
+  input: PolygonData,
+  callback: MapCoordsCallback,
+): PolygonData;
+export function mapCoords(
+  input: MultiPointData,
+  callback: MapCoordsCallback,
+): MultiPointData;
+export function mapCoords(
+  input: MultiLineStringData,
+  callback: MapCoordsCallback,
+): MultiLineStringData;
+export function mapCoords(
+  input: MultiPolygonData,
+  callback: MapCoordsCallback,
+): MultiPolygonData;
+
+// TODO: ideally I could use <T extends GeoArrowType> here...
+export function mapCoords(
+  input: GeoArrowData,
+  callback: MapCoordsCallback,
+): GeoArrowData {
+  if (isPointData(input)) {
+    return mapCoords0(input, callback);
+  }
+  if (isLineStringData(input)) {
+    return mapCoords1(input, callback);
+  }
+  if (isPolygonData(input)) {
+    return mapCoords2(input, callback);
+  }
+  if (isMultiPointData(input)) {
+    return mapCoords1(input, callback);
+  }
+  if (isMultiLineStringData(input)) {
+    return mapCoords2(input, callback);
+  }
+  if (isMultiPolygonData(input)) {
+    return mapCoords3(input, callback);
+  }
+
+  assertFalse();
+}
+
+export function mapCoords0<T extends Point>(
+  input: arrow.Data<T>,
+  callback: MapCoordsCallback,
+): arrow.Data<T> {
+  assert(input.type.listSize === 2, "expected 2D");
+  const coordsData = getPointChild(input);
+  const flatCoords = coordsData.values;
+
+  const outputCoords = new Float64Array(flatCoords.length);
+  for (let coordIdx = 0; coordIdx < input.length; coordIdx++) {
+    const x = flatCoords[coordIdx * 2];
+    const y = flatCoords[coordIdx * 2 + 1];
+    const [newX, newY] = callback(x, y);
+    outputCoords[coordIdx * 2] = newX;
+    outputCoords[coordIdx * 2 + 1] = newY;
+  }
+
+  const newCoordsData = arrow.makeData({
+    type: coordsData.type,
+    length: coordsData.length,
+    nullCount: coordsData.nullCount,
+    nullBitmap: coordsData.nullBitmap,
+    data: outputCoords,
+  });
+
+  return arrow.makeData({
+    type: input.type,
+    length: input.length,
+    nullCount: input.nullCount,
+    nullBitmap: input.nullBitmap,
+    child: newCoordsData,
+  });
+}
+
+/**
+ * NOTE: the callback must be infallible as this does not take geometry validity
+ * into effect for operating on coords
+ */
+export function mapCoords1<T extends LineString | MultiPoint>(
+  input: arrow.Data<T>,
+  callback: MapCoordsCallback,
+): arrow.Data<T> {
+  const pointData = getLineStringChild(input);
+  const newPointData = mapCoords0(pointData, callback);
+
+  return arrow.makeData({
+    type: input.type,
+    length: input.length,
+    nullCount: input.nullCount,
+    nullBitmap: input.nullBitmap,
+    child: newPointData,
+    valueOffsets: input.valueOffsets,
+  });
+}
+
+/**
+ * NOTE: the callback must be infallible as this does not take geometry validity
+ * into effect for operating on coords
+ */
+export function mapCoords2<T extends Polygon | MultiLineString>(
+  input: arrow.Data<T>,
+  callback: MapCoordsCallback,
+): arrow.Data<T> {
+  const linestringData = getPolygonChild(input);
+  const newLinestringData = mapCoords1(linestringData, callback);
+
+  return arrow.makeData({
+    type: input.type,
+    length: input.length,
+    nullCount: input.nullCount,
+    nullBitmap: input.nullBitmap,
+    child: newLinestringData,
+    valueOffsets: input.valueOffsets,
+  });
+}
+
+/**
+ * NOTE: the callback must be infallible as this does not take geometry validity
+ * into effect for operating on coords
+ */
+export function mapCoords3<T extends MultiPolygon>(
+  input: arrow.Data<T>,
+  callback: MapCoordsCallback,
+): arrow.Data<T> {
+  const polygonData = getMultiPolygonChild(input);
+  const newPolygonData = mapCoords2(polygonData, callback);
+
+  return arrow.makeData({
+    type: input.type,
+    length: input.length,
+    nullCount: input.nullCount,
+    nullBitmap: input.nullBitmap,
+    child: newPolygonData,
+    valueOffsets: input.valueOffsets,
+  });
+}

src/algorithm/index.ts

@@ -1,3 +1,6 @@
 export { area, signedArea } from "./area.js";
 export { earcut } from "./earcut.js";
+export { mapCoords } from "./coords.js";
+export { reproject } from "./proj.js";
+export { totalBounds } from "./total-bounds.js";
 export { windingDirection, modifyWindingDirection } from "./winding.js";

src/algorithm/proj.ts

@@ -0,0 +1,53 @@
+import * as arrow from "apache-arrow";
+import proj4 from "proj4";
+import { GeoArrowType } from "../type";
+import { mapCoords } from "./coords";
+
+/**
+ * Reproject using proj4
+ */
+export function reproject<T extends GeoArrowType>(
+  input: arrow.Data<T>,
+  fromProjection: string,
+  toProjection: string,
+): arrow.Data<T>;
+export function reproject<T extends GeoArrowType>(
+  input: arrow.Vector<T>,
+  fromProjection: string,
+  toProjection: string,
+): arrow.Vector<T>;
+
+export function reproject<T extends GeoArrowType>(
+  input: arrow.Data<T> | arrow.Vector<T>,
+  fromProjection: string,
+  toProjection: string,
+): arrow.Data<T> | arrow.Vector<T> {
+  const projectionFn = proj4(fromProjection, toProjection);
+  if ("data" in input) {
+    return new arrow.Vector(
+      input.data.map((data) => reprojectData(data, projectionFn)),
+    );
+  }
+
+  return reprojectData(input, projectionFn);
+}
+
+/**
+ * Reproject a single Data instance
+ */
+function reprojectData<T extends GeoArrowType>(
+  input: arrow.Data<T>,
+  projectionFn: proj4.Converter,
+): arrow.Data<T> {
+  // Avoid extra object creation
+  const stack = [0, 0];
+  const callback = (x: number, y: number) => {
+    stack[0] = x;
+    stack[1] = y;
+    return projectionFn.forward(stack) as [number, number];
+  };
+
+  // @ts-expect-error I have a mismatch between generic T extends GeoArrowType
+  // and concrete GeoArrowData typing
+  return mapCoords(input, callback);
+}

src/algorithm/total-bounds.ts

@@ -99,34 +99,16 @@ function totalBoundsNest0(vector: PointVector): Bbox {
 }
 
 function totalBoundsNest1(vector: LineStringVector): Bbox {
-  const bbox = new Bbox();
-  for (const data of vector.data) {
-    const pointData = getLineStringChild(data);
-    bbox.updateBbox(coordsBbox(pointData));
-  }
-
-  return bbox;
+  const pointVector = getLineStringChild(vector);
+  return totalBoundsNest0(pointVector);
 }
 
 function totalBoundsNest2(vector: PolygonVector): Bbox {
-  const bbox = new Bbox();
-  for (const data of vector.data) {
-    const lineStringData = getPolygonChild(data);
-    const pointData = getLineStringChild(lineStringData);
-    bbox.updateBbox(coordsBbox(pointData));
-  }
-
-  return bbox;
+  const lineStringVector = getPolygonChild(vector);
+  return totalBoundsNest1(lineStringVector);
 }
 
 function totalBoundsNest3(vector: MultiPolygonVector): Bbox {
-  const bbox = new Bbox();
-  for (const data of vector.data) {
-    const polygonData = getMultiPolygonChild(data);
-    const lineStringData = getPolygonChild(polygonData);
-    const pointData = getLineStringChild(lineStringData);
-    bbox.updateBbox(coordsBbox(pointData));
-  }
-
-  return bbox;
+  const polygonVector = getMultiPolygonChild(vector);
+  return totalBoundsNest2(polygonVector);
 }

src/algorithm/utils/assert.ts

@@ -0,0 +1,9 @@
+export function assert(condition: boolean, message?: string) {
+  if (!condition) {
+    throw new Error(`assertion failed ${message}`);
+  }
+}
+
+export function assertFalse(): never {
+  throw new Error(`assertion failed`);
+}

src/data.ts

@@ -20,6 +20,13 @@ export type PolygonData = arrow.Data<Polygon>;
 export type MultiPointData = arrow.Data<MultiPoint>;
 export type MultiLineStringData = arrow.Data<MultiLineString>;
 export type MultiPolygonData = arrow.Data<MultiPolygon>;
+export type GeoArrowData =
+  | PointData
+  | LineStringData
+  | PolygonData
+  | MultiPointData
+  | MultiLineStringData
+  | MultiPolygonData;
 
 export function isPointData(data: arrow.Data): data is PointData {
   return isPoint(data.type);

src/type.ts

@@ -1,9 +1,12 @@
 import * as arrow from "apache-arrow";
 
-export type InterleavedCoord = arrow.FixedSizeList<arrow.Float64>;
+// Note: this apparently has to be arrow.Float and not arrow.Float64 to ensure
+// that recreating a data instance with arrow.makeData type checks using the
+// input's data type.
+export type InterleavedCoord = arrow.FixedSizeList<arrow.Float>;
 export type SeparatedCoord = arrow.Struct<{
-  x: arrow.Float64;
-  y: arrow.Float64;
+  x: arrow.Float;
+  y: arrow.Float;
 }>;
 // TODO: support separated coords
 export type Coord = InterleavedCoord; // | SeparatedCoord;
@@ -13,6 +16,13 @@ export type Polygon = arrow.List<arrow.List<Coord>>;
 export type MultiPoint = arrow.List<Coord>;
 export type MultiLineString = arrow.List<arrow.List<Coord>>;
 export type MultiPolygon = arrow.List<arrow.List<arrow.List<Coord>>>;
+export type GeoArrowType =
+  | Point
+  | LineString
+  | Polygon
+  | MultiPoint
+  | MultiLineString
+  | MultiPolygon;
 
 /** Check that the given type is a Point data type */
 export function isPoint(type: arrow.DataType): type is Point {

src/vector.ts

@@ -20,6 +20,13 @@ export type PolygonVector = arrow.Vector<Polygon>;
 export type MultiPointVector = arrow.Vector<MultiPoint>;
 export type MultiLineStringVector = arrow.Vector<MultiLineString>;
 export type MultiPolygonVector = arrow.Vector<MultiPolygon>;
+export type GeoArrowVector =
+  | PointVector
+  | LineStringVector
+  | PolygonVector
+  | MultiPointVector
+  | MultiLineStringVector
+  | MultiPolygonVector;
 
 export function isPointVector(vector: arrow.Vector): vector is PointVector {
   return isPoint(vector.type);

tests/algorithm/proj.test.ts

@@ -0,0 +1,23 @@
+import { describe, expect, it } from "vitest";
+import { testPointData } from "../util/point";
+import { reproject } from "../../src/algorithm";
+
+import proj4 from "proj4";
+
+describe("reproject", (t) => {
+  it("should reproject point array", () => {
+    const pointData = testPointData();
+    const reprojected = reproject(pointData, "EPSG:4326", "EPSG:3857");
+
+    const expected1 = proj4("EPSG:4326", "EPSG:3857", [1, 2]);
+    const expected2 = proj4("EPSG:4326", "EPSG:3857", [3, 4]);
+    const expected3 = proj4("EPSG:4326", "EPSG:3857", [5, 6]);
+
+    expect(reprojected.children[0].values[0]).toBeCloseTo(expected1[0]);
+    expect(reprojected.children[0].values[1]).toBeCloseTo(expected1[1]);
+    expect(reprojected.children[0].values[2]).toBeCloseTo(expected2[0]);
+    expect(reprojected.children[0].values[3]).toBeCloseTo(expected2[1]);
+    expect(reprojected.children[0].values[4]).toBeCloseTo(expected3[0]);
+    expect(reprojected.children[0].values[5]).toBeCloseTo(expected3[1]);
+  });
+});