|
| 1 | +/* |
| 2 | + * Complex logarithm projection |
| 3 | + * |
| 4 | + * Based on the following papers by Joachim Böttger et al.: |
| 5 | + * - Detail‐In‐Context Visualization for Satellite Imagery (2008) (https://doi.org/10.1111/j.1467-8659.2008.01156.x) |
| 6 | + * - Complex Logarithmic Views for Small Details in Large Contexts (2006) (https://doi.org/10.1109/TVCG.2006.126) |
| 7 | + * |
| 8 | + * Implemented for d3 by Matthias Albrecht and Jochen Görtler (2019) |
| 9 | + * |
| 10 | + */ |
| 11 | + |
| 12 | +import { geoProjectionMutator as projectionMutator, geoAzimuthalEqualAreaRaw as azimuthalEqualAreaRaw } from "d3-geo"; |
| 13 | +import { abs, sin, cos, pi, exp, atan2 } from "./math.js"; |
| 14 | +import { complexMul, complexLogHypot } from "./complex.js"; |
| 15 | +import { default as clipPolygon } from "./clip/polygon.js"; |
| 16 | + |
| 17 | +// Default planar projection and cutoff latitude, see below for an explanation of these settings. |
| 18 | +var DEFAULT_PLANAR_PROJECTION_RAW = azimuthalEqualAreaRaw; |
| 19 | +var DEFAULT_CUTOFF_LATITUDE = -0.05; |
| 20 | + |
| 21 | +// Offset used to prevent logarithm of 0. |
| 22 | +var CARTESIAN_OFFSET = 1e-10; |
| 23 | + |
| 24 | +// Projection parameters for the default 960x500 projection area. |
| 25 | +var DEFAULT_PROJECTION_PARAMS = { |
| 26 | + angle: 90, |
| 27 | + center: [0, 5.022570623227068], |
| 28 | + scale: 79.92959180396787, |
| 29 | + translate: [479.9999905630355, 250.35977064160338] |
| 30 | +} |
| 31 | + |
| 32 | +// Vertices of the clipping polygon in spherical coordinates. |
| 33 | +// It contains the whole world except a small strip along longitude 0/180 crossing the south pole. |
| 34 | +var CLIP_POLY_SPHERICAL = [ |
| 35 | + [-180, -1e-4], |
| 36 | + [180, -1e-4], |
| 37 | + [1e-4, DEFAULT_CUTOFF_LATITUDE], |
| 38 | + [-1e-4, DEFAULT_CUTOFF_LATITUDE] |
| 39 | +] |
| 40 | + |
| 41 | +// Clipping polygon precision. |
| 42 | +var N_SIDE = 5; |
| 43 | +var N_BOTTOM = 50; |
| 44 | + |
| 45 | + |
| 46 | +export function complexLogRaw(planarProjectionRaw = DEFAULT_PLANAR_PROJECTION_RAW) { |
| 47 | + function forward(lambda, phi) { |
| 48 | + // Project on plane. |
| 49 | + // Interpret projected point on complex plane. |
| 50 | + var aziComp = planarProjectionRaw(lambda, phi); |
| 51 | + |
| 52 | + // Rotate by -90 degrees in complex plane so the following complex log projection will be horizontally centered |
| 53 | + aziComp = complexMul(aziComp, [cos(-pi / 2), sin(-pi / 2)]); |
| 54 | + |
| 55 | + // Small offset to prevent logarithm of 0. |
| 56 | + if (aziComp[0] == 0 && aziComp[1] == 0) { |
| 57 | + aziComp[0] += CARTESIAN_OFFSET; |
| 58 | + aziComp[1] += CARTESIAN_OFFSET; |
| 59 | + } |
| 60 | + |
| 61 | + // Apply complex logarithm. |
| 62 | + var logComp = [complexLogHypot(aziComp[0], aziComp[1]), atan2(aziComp[1], aziComp[0])]; |
| 63 | + |
| 64 | + return logComp; |
| 65 | + } |
| 66 | + |
| 67 | + function invert(x, y) { |
| 68 | + // Inverse complex logarithm (complex exponential function). |
| 69 | + var invLogComp = [exp(x) * cos(y), exp(x) * sin(y)]; |
| 70 | + |
| 71 | + // Undo rotation. |
| 72 | + invLogComp = complexMul(invLogComp, [cos(pi / 2), sin(pi / 2)]); |
| 73 | + |
| 74 | + // Invert azimuthal equal area. |
| 75 | + return planarProjectionRaw.invert(invLogComp[0], invLogComp[1]); |
| 76 | + } |
| 77 | + |
| 78 | + forward.invert = invert; |
| 79 | + return forward; |
| 80 | +} |
| 81 | + |
| 82 | + |
| 83 | +export default function(planarProjectionRaw = DEFAULT_PLANAR_PROJECTION_RAW, cutoffLatitude = DEFAULT_CUTOFF_LATITUDE) { |
| 84 | + var mutator = projectionMutator(complexLogRaw); |
| 85 | + var projection = mutator(planarProjectionRaw); |
| 86 | + |
| 87 | + // Projection used to project onto the complex plane. |
| 88 | + projection.planarProjectionRaw = function(_) { |
| 89 | + return arguments.length ? clipped(mutator(planarProjectionRaw = _)) : planarProjectionRaw; |
| 90 | + } |
| 91 | + |
| 92 | + // Latitude relative to the projection center at which to cutoff/clip the projection, lower values result in more detail around the projection center. |
| 93 | + // Value must be < 0 because complex log projects the origin to infinity. |
| 94 | + projection.cutoffLatitude = function(_) { |
| 95 | + return arguments.length ? (cutoffLatitude = _, clipped(mutator(planarProjectionRaw))) : cutoffLatitude; |
| 96 | + } |
| 97 | + |
| 98 | + function clipped(projection) { |
| 99 | + var angle = projection.angle(); |
| 100 | + var scale = projection.scale(); |
| 101 | + var center = projection.center(); |
| 102 | + var translate = projection.translate(); |
| 103 | + var rotate = projection.rotate(); |
| 104 | + |
| 105 | + projection |
| 106 | + .angle(DEFAULT_PROJECTION_PARAMS.angle) |
| 107 | + .scale(1) |
| 108 | + .center([0, 0]) |
| 109 | + .rotate([0, 0]) |
| 110 | + .translate([0, 0]) |
| 111 | + .preclip(); |
| 112 | + |
| 113 | + // These are corner vertices of a rectangle in the projected complex log view. |
| 114 | + var topLeft = projection(CLIP_POLY_SPHERICAL[0]); |
| 115 | + var topRight = projection(CLIP_POLY_SPHERICAL[1]); |
| 116 | + var bottomRight = projection([CLIP_POLY_SPHERICAL[2][0], cutoffLatitude]); |
| 117 | + var bottomLeft = projection([CLIP_POLY_SPHERICAL[3][0], cutoffLatitude]); |
| 118 | + var width = abs(topRight[0] - topLeft[0]); |
| 119 | + var height = abs(bottomRight[1] - topRight[1]); |
| 120 | + |
| 121 | + // Prevent overlapping polygons that result from paths that go from one side to the other, |
| 122 | + // so cut along 180°/-180° degree line (left and right in complex log projected view). |
| 123 | + // This means cutting against a rectangular shaped polygon in the projected view. |
| 124 | + // The following generator produces a polygon that is shaped like this: |
| 125 | + // |
| 126 | + // Winding order: ==> |
| 127 | + // |
| 128 | + // ******************| |
| 129 | + // | | |
| 130 | + // | | |
| 131 | + // | | |
| 132 | + // | | |
| 133 | + // | | |
| 134 | + // |------------------ |
| 135 | + // |
| 136 | + // N_SIDE determines how many vertices to insert along the sides (marked as | above). |
| 137 | + // N_BOTTOM determines how many vertices to insert along the bottom (marked as - above). |
| 138 | + // |
| 139 | + // The resulting polygon vertices are back-projected to spherical coordinates. |
| 140 | + var polygon = { |
| 141 | + type: "Polygon", |
| 142 | + coordinates: [ |
| 143 | + [ |
| 144 | + topLeft, |
| 145 | + ...Array.from({length: N_SIDE}, (_, t) => [bottomRight[0], bottomRight[1] - height * (N_SIDE- t) / N_SIDE]), |
| 146 | + ...Array.from({length: N_BOTTOM}, (_, t) => [bottomRight[0] - width * t / N_BOTTOM, bottomRight[1]]), |
| 147 | + ...Array.from({length: N_SIDE}, (_, t) => [bottomLeft[0], bottomLeft[1] - height * t / N_SIDE]), |
| 148 | + topLeft |
| 149 | + ].map(point => projection.invert(point)) |
| 150 | + ] |
| 151 | + }; |
| 152 | + |
| 153 | + return projection |
| 154 | + .angle(angle) |
| 155 | + .scale(scale) |
| 156 | + .center(center) |
| 157 | + .translate(translate) |
| 158 | + .rotate(rotate) |
| 159 | + .preclip(clipPolygon(polygon)); |
| 160 | + } |
| 161 | + |
| 162 | + // The following values are for the default 960x500 projection area |
| 163 | + return clipped(projection) |
| 164 | + .angle(DEFAULT_PROJECTION_PARAMS.angle) |
| 165 | + .center(DEFAULT_PROJECTION_PARAMS.center) |
| 166 | + .scale(DEFAULT_PROJECTION_PARAMS.scale) |
| 167 | + .translate(DEFAULT_PROJECTION_PARAMS.translate); |
| 168 | +} |
0 commit comments