orbis-core/src/projection.ts

import * as d3geo from 'd3-geo';
import * as d3geoProjection from 'd3-geo-projection';
import * as d3geoPolygon from 'd3-geo-polygon';
import { PNG } from 'pngjs';
import * as Png from './png';

export type Point = [number, number]
export type Bounds = [Point, Point]

type Projection = [string, ...unknown[]]

type Rect = {
  width: number,
  height: number,
}

export type ProjectOptions = {
  bounds: Bounds,
  wrapAround: boolean,
  outputSize?: Partial<Rect>,
  outputPadding?: {
    x?: number,
    y?: number,
  }
}

const geoProjectionNamePrefix = 'geo';

type GeoProjectionFactory = any

const getProjectionFunction = (projectionFunctionName: string): GeoProjectionFactory => {
  if (projectionFunctionName in d3geoPolygon) {
    return (d3geoPolygon as Record<string, unknown>)[projectionFunctionName];
  }
  if (projectionFunctionName in d3geoProjection) {
    return (d3geoProjection as Record<string, unknown>)[projectionFunctionName];
  }
  if (projectionFunctionName in d3geo) {
    return (d3geo as Record<string, unknown>)[projectionFunctionName];
  }

  const properProjectionName = projectionFunctionName.slice(
    projectionFunctionName.indexOf(geoProjectionNamePrefix)
    + geoProjectionNamePrefix.length,
  );

  throw new Error(`Unknown projection: ${properProjectionName}`);
};

const buildProjection = (projection: Projection) => {
  const [projectionName, ...projectionArgs] = projection;
  // eslint-disable-next-line @typescript-eslint/no-unsafe-assignment
  const projectionFunction = getProjectionFunction(`${geoProjectionNamePrefix}${projectionName}`);
  // eslint-disable-next-line @typescript-eslint/no-unsafe-call
  return projectionFunction(projectionArgs) as d3geo.GeoProjection;
};

const getCenter = (bounds: Bounds): Point => {
  const [nw, se] = bounds;
  const [nwLng, nwLat] = nw;
  const [seLng, seLat] = se;
  let seLngNorm = seLng;

  while (seLngNorm < nwLng) {
    seLngNorm += 360;
  }

  return [
    ((((nwLng + 180) + (seLngNorm + 180)) / 2) % 360) - 180,
    (nwLat + seLat) / 2,
  ];
};

const transformGeoProjection = (
  geoProjection: d3geo.GeoProjection,
  options: ProjectOptions,
  inputWidth: number,
) => {
  const center = getCenter(options.bounds);
  const transformedGeoProjection = geoProjection
    .reflectX(false)
    .reflectY(false)
    .rotate([-center[0], -center[1]]);

  const boundsGeoJson: GeoJSON.Polygon = {
    type: 'Polygon',
    coordinates: [
      [
        [options.bounds[0][0], options.bounds[0][1]],
        [options.bounds[1][0], options.bounds[0][1]],
        [options.bounds[1][0], options.bounds[1][1]],
        [options.bounds[0][0], options.bounds[1][1]],
      ],
    ],
  };

  if (typeof options.outputSize?.width === 'number') {
    if (typeof options.outputSize.height === 'number') {
      const paddingX = options.outputPadding?.x ?? 0;
      const paddingY = options.outputPadding?.y ?? 0;
      return transformedGeoProjection.fitExtent(
        [
          [paddingX, paddingY],
          [options.outputSize.width - paddingX, options.outputSize.height - paddingY],
        ],
        boundsGeoJson,
      );
    }
    return transformedGeoProjection.fitWidth(
      options.outputSize.width,
      boundsGeoJson,
    );
  }
  if (typeof options.outputSize?.height === 'number') {
    return transformedGeoProjection.fitHeight(
      options.outputSize.height,
      boundsGeoJson,
    );
  }

  return transformedGeoProjection.fitWidth(inputWidth, boundsGeoJson);
};

const computeLatitudeSensitiveProjectionVerticalBounds = (
  transformedGeoProjection: d3geo.GeoProjection,
  bounds: Bounds,
  outputWidthRaw: number,
) => {
  // web mercator clipping
  // const maxAbsLatitude
  //   = ((2 * Math.atan(Math.E ** Math.PI) - (Math.PI / 2)) / (Math.PI * 2)) * 360;
  const maxAbsLatitude = 85.0511287798066;
  const adjustedNwBoundsRaw = transformedGeoProjection([
    bounds[0][0],
    Math.min(bounds[0][1], maxAbsLatitude),
  ]);
  const adjustedSeBoundsRaw = transformedGeoProjection([
    bounds[1][0],
    Math.max(bounds[1][1], -maxAbsLatitude),
  ]);

  if (!(adjustedNwBoundsRaw && adjustedSeBoundsRaw)) {
    return null;
  }

  return [
    Math.round(outputWidthRaw),
    Math.round(adjustedSeBoundsRaw[1] - adjustedNwBoundsRaw[1]),
  ];
};

const computeLongitudeSensitiveProjectionHorizontalBounds = (
  transformedGeoProjection: d3geo.GeoProjection,
  bounds: Bounds,
  outputHeightRaw: number,
) => {
  const adjustedWBoundsRaw = transformedGeoProjection([
    bounds[0][0],
    (bounds[0][1] + bounds[1][1]) / 2,
  ]);
  const adjustedEBoundsRaw = transformedGeoProjection([
    bounds[1][0],
    (bounds[0][1] + bounds[1][1]) / 2,
  ]);

  if (!(adjustedWBoundsRaw && adjustedEBoundsRaw)) {
    return null;
  }

  return [
    Math.round(adjustedEBoundsRaw[0] - adjustedWBoundsRaw[0]),
    Math.round(outputHeightRaw),
  ];
};

const computeOutputBounds = (
  transformedGeoProjection: d3geo.GeoProjection,
  projectionName: string,
  bounds: Bounds,
  outputSize?: Partial<Rect>,
) => {
  // TODO: what would it be for polygonal projections?

  if (
    typeof outputSize?.width === 'number'
    && typeof outputSize.height === 'number'
  ) {
    return [
      outputSize.width,
      outputSize.height,
    ];
  }

  const nwBoundsRaw = transformedGeoProjection(bounds[0]);
  const seBoundsRaw = transformedGeoProjection(bounds[1]);

  if (!(nwBoundsRaw && seBoundsRaw)) {
    return null;
  }

  const outputWidthRaw = seBoundsRaw[0] - nwBoundsRaw[0];
  const outputHeightRaw = seBoundsRaw[1] - nwBoundsRaw[1];

  switch (projectionName) {
  case 'Mercator':
    return computeLatitudeSensitiveProjectionVerticalBounds(
      transformedGeoProjection,
      bounds,
      outputWidthRaw,
    );
  case 'Robinson':
  case 'Sinusoidal':
    return computeLongitudeSensitiveProjectionHorizontalBounds(
      transformedGeoProjection,
      bounds,
      outputHeightRaw,
    );
  default:
    break;
  }

  return [
    Math.round(outputWidthRaw),
    Math.round(outputHeightRaw),
  ];
};

export const project = async (
  pngInput: string | Buffer | PNG,
  projection: Projection,
  options = {
    wrapAround: false,
    bounds: [[-180, 90], [180, -90]],
    outputSize: undefined,
  } as ProjectOptions,
) => {
  const inputPng = await Png.load(pngInput);
  const geoProjection = buildProjection(projection);
  const [projectionName] = projection;

  const transformedGeoProjection = transformGeoProjection(geoProjection, options, inputPng.width);
  if (!transformedGeoProjection.invert) {
    throw new Error(`No invert() function for projection "${projectionName}"`);
  }

  const outputBounds = computeOutputBounds(
    transformedGeoProjection,
    projectionName,
    options.bounds,
    options.outputSize,
  );

  if (!outputBounds) {
    throw new Error(
      `Cannot compute bounds, possibly unimplemented. Check logic of "${projectionName}" projection.`,
    );
  }

  const [outputWidth, outputHeight] = outputBounds;
  // outputWidth = options.outputSize?.width ?? 461;
  // outputHeight = options.outputSize?.height ?? 461;

  const { data: inputData } = inputPng;
  const outputPng = new PNG({
    width: outputWidth,
    height: outputHeight,
  });
  const { data: outputData } = outputPng;

  let i = 0;
  for (let y = 0; y < outputHeight; y += 1) {
    for (let x = 0; x < outputWidth; x += 1) {
      const projected = transformedGeoProjection.invert([x, y]);
      if (projected) {
        const reversed = transformedGeoProjection(projected);
        if (reversed) {
          const isSamePoint = (
            Math.abs(reversed[0] - x) < 0.5
            && Math.abs(reversed[1] - y) < 0.5
          );

          if (isSamePoint) {
            const [lambda, phi] = projected;
            if (!(lambda > 180 || lambda < -180 || phi > 90 || phi < -90)) {
              // eslint-disable-next-line no-bitwise
              const q = (((90 - phi) / 180) * inputPng.height | 0) * inputPng.width
                // eslint-disable-next-line no-bitwise
                + (((180 + lambda) / 360) * inputPng.width | 0) << 2;
              outputData[i] = inputData[q];
              outputData[i + 1] = inputData[q + 1];
              outputData[i + 2] = inputData[q + 2];
              outputData[i + 3] = inputData[q + 3];
            }
          }
        }
      }
      i += 4;
    }
  }

  return PNG.sync.write(outputPng);
};