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';
import { ProjectBounds } from './common';
import { getCountryGeometry } from './geo';

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

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

export type ProjectOptions = {
  bounds: ProjectBounds,
  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 loadProjection = (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 getMidpoint = (a: GeoJSON.Position, b: GeoJSON.Position) => [
  (a[0] + b[0]) / 2,
  (a[1] + b[1]) / 2,
];

const getCentroid = (boundsRaw: d3geo.GeoGeometryObjects): GeoJSON.Position => {
  if (boundsRaw.type === 'Polygon') {
    return boundsRaw.coordinates[0].reduce(getMidpoint);
  }

  if (boundsRaw.type === 'MultiPolygon') {
    return boundsRaw.coordinates
      .map((polygon) => polygon[0].reduce(getMidpoint))
      .reduce(getMidpoint);
  }

  throw new Error(`Invalid type: ${boundsRaw.type}`);
};

const buildGeoProjection = (
  projection: Projection,
  boundsGeoJson: d3geo.GeoGeometryObjects,
  options: ProjectOptions,
  inputWidth: number,
) => {
  const geoProjection = loadProjection(projection);
  const center2 = d3geo.geoCentroid(boundsGeoJson); // why is d3geo's centroid different than our simple centroid function?
  const center = getCentroid(boundsGeoJson);

  console.log(center, center2);
  const transformedGeoProjection = geoProjection
    .reflectX(false)
    .reflectY(false)
    .rotate([-center[0], -center[1]]);

  // fixme: fitX() methods don't modify scale, only the dimensions?
  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 computeOutputBounds = (
  transformedGeoProjection: d3geo.GeoProjection,
  _projectionName: string,
  bounds: d3geo.GeoGeometryObjects,
) => {
  const path = d3geo.geoPath(transformedGeoProjection);

  const bbox = path.bounds(bounds);
  // console.log(bbox);

  return [
    Math.round(Math.abs(bbox[1][0] - bbox[0][0])),
    Math.round(Math.abs(bbox[1][1] - bbox[0][1])),
  ];
};

// const computeOutputBounds = (
//   transformedGeoProjection: d3geo.GeoProjection,
//   projectionName: string,
//   bounds: d3geo.GeoGeometryObjects,
//   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 path = d3.geoPath(transformedGeoProjection);
//
//   path.bounds(bounds)
//   // TODO better way to compute width and height directly from d3?
//
//   let extremeProjectedBoundsRaw: Bounds | null | undefined;
//   if (bounds.type === 'Polygon') {
//     extremeProjectedBoundsRaw = bounds.coordinates[0].reduce(
//       (theBounds, point) => {
//         const projected = transformedGeoProjection(point as [number, number]);
//
//         if (projected === null) {
//           return null;
//         }
//
//         if (theBounds === null) {
//           return [
//             projected,
//             projected,
//           ];
//         }
//
//         return [
//           [
//             Math.min(theBounds[0][0], projected[0]),
//             Math.max(theBounds[0][1], projected[1]),
//           ],
//           [
//             Math.max(theBounds[1][0], projected[0]),
//             Math.min(theBounds[1][1], projected[1]),
//           ],
//         ];
//       },
//       null as unknown as (GeoJSON.Position[] | null),
//     );
//   }
//
//   if (!extremeProjectedBoundsRaw) {
//     return null;
//   }
//
//   const [nwBoundsRaw, seBoundsRaw] = extremeProjectedBoundsRaw;
//   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),
//   ];
// };

const determineGeoJsonBounds = async (bounds: ProjectBounds): Promise<d3geo.GeoGeometryObjects> => {
  if (bounds.type === 'country') {
    const countryGeometry = await getCountryGeometry(
      bounds.value,
      { mergeBoundingBoxes: true },
    );
    return countryGeometry.geometry;
  }

  if (bounds.type === 'geojson') {
    return bounds.value;
  }

  throw new Error(`Invalid bounds type: ${(bounds as Record<string, string>).type}`);
};

export const project = async (
  pngInput: string | Buffer | PNG,
  projection: Projection,
  options = {
    wrapAround: false,
    bounds: {
      type: 'geojson',
      value: {
        type: 'Sphere',
      },
    },
    outputSize: undefined,
  } as ProjectOptions,
) => {
  console.log(options);
  const inputPng = await Png.load(pngInput);
  const theBounds = await determineGeoJsonBounds(options.bounds);

  const transformedGeoProjection = buildGeoProjection(
    projection,
    theBounds,
    options,
    inputPng.width,
  );

  const [projectionName] = projection;
  if (!transformedGeoProjection.invert) {
    throw new Error(`No invert() function for projection "${projectionName}"`);
  }

  const outputBounds = computeOutputBounds(
    transformedGeoProjection,
    projectionName,
    theBounds,
  );

  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) {
        if (options.wrapAround) {
          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];
          }
        } else {
          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);
};