30:T4c9,md`# Tissot's indicatrix When thinking about [Tissot's indicatrix](https://en.wikipedia.org/wiki/Tissot%27s_indicatrix), probably most people imagine it as a tool to visualize spatial change in distortion, illustrated by more or less elliptical shapes placed across a map. In most instances, however, these visualizations [don't show true indicatrices](https://observablehq.com/@d3/tissots-indicatrix), but projected circles placed around a reference point at a given radius, that give an intuitive representation of areal and shape distortion on a map. So what are *true Tissot's indicatrices* then? [Nicolas Auguste Tissot](https://en.wikipedia.org/wiki/Nicolas_Auguste_Tissot) published his classic analysis on the distortion on maps in 1859 and 1881. The basic idea is that the intersection of any two lines on the Earth is represented on the flat map with an intersection at the same or a different angle. He proved that at almost every point on the Earth, there's a right angle intersection of two lines in *some* direction which are also shown at right angles on the map. All the other intersections at that point will not intersect at the same angle on the map, unless the map is conformal, at least at that point.`31:T96c,{ const w = Math.min(640, width), height = Math.min(320, width / 2); const svg = d3.select(DOM.svg(w, height)); const map = svg.append("g") .attr("fill", "#ececec") .attr("stroke", "#dcdcdc") const plot = svg.append("g") .attr("transform", `translate(${w / 2}, ${height / 2}) scale(100, 100)`) .attr("fill", "none") .attr("stroke", "black") .attr("stroke-width", 1 / 100); const dasharray = `${5 / 100} ${2 / 100}`; // Note: // Only use projections with straight parallels. // Axis rotation is not taken into account in this example.. const projection = d3.geoSinusoidal() .translate([w - 80, 60]) .scale(25); const path = d3.geoPath().projection(projection); const xy = projection([lon, lat]); // background plot .append("circle") .attr("r", 1) .attr("fill", "none") .attr("stroke", "grey") .attr("stroke-dasharray", dasharray); map.append("path") .attr("d", path(land)) .attr("stroke", "none"); map.append("path") .attr("d", path(sphere)) .attr("fill", "none"); map.append("circle") .attr("r", 2) .attr("cx", xy[0]) .attr("cy", xy[1]) .attr("fill", "#cc0000"); // Indicatrix const tissot = geoTissot().projection(projection); const latClamped = Math.max(-90 + 1e-3, Math.min(90 - 1e-3, lat)); const indicatrix = tissot({ type: "Point", coordinates: [lon, latClamped] })[0]; if (!indicatrix) return svg.node(); let { s, h, k, a, b, theta } = indicatrix; let rotation = Math.asin(s / (h * k)) * Math.sign(theta) * deg || 90; if (lon === 0 || lat === 0) theta = Math.PI / 4; // parallel phi plot .append("line") .attr("x1", -k) .attr("x2", k); // meridian lambda plot .append("line") .attr("transform", `rotate(${rotation})`) .attr("x1", -h) .attr("x2", h); // indicatrix const ellipse = plot.append("g") .attr("transform", `rotate(${theta * deg})`); ellipse .append("line") .attr("x1", -a) .attr("x2", a) .attr("stroke-dasharray", dasharray); ellipse .append("line") .attr("y1", -b) .attr("y2", b) .attr("stroke-dasharray", dasharray); ellipse .append("ellipse") .attr("rx", a) .attr("ry", b); return html`

Tissot's indicatrix for a sinusoidal projection

${svg.node()} `; }32:T491,md`Tissot showed this relationship graphically with a special *ellipse of distortion* called an *indicatrix* by computing partial derivatives at specific points. An infinitely small circle on the Earth projects as an infinitely small, but perfect, ellipse on any map projection. If the projection is conformal, the ellipse is a circle, an ellipse of zero eccentricity. Otherwise, the ellipse has a major axis ${tex`a`} and minor axis ${tex`b`} which are directly related to the *maximum angular deformation* ${tex`\omega`} and to the *areal scale distortion* ${tex`s`}. The image below shows *areal scale distortion* on a map with [Mercator projection](https://en.wikipedia.org/wiki/Mercator_projection). It preserves angles around all locations, so all indicatrices are shown as circles. However, the *areal scale distortion* increases as distance away from the line of true scale increases (the equator on an unrotated projection). For example, an area located 60° N or S form the line of true scale will appear four times larger on the map than an area of the same size located at the line of true scale. (Make sure to drag the map to explore the effect)
`33:Tb08,{ const size = Math.min(600, width - 40); const context = DOM.context2d(size + 40, size); const mercator = d3.geoMercator() //.rotate([110, -12, 20]) .fitSize([size, size], sphere); const path = d3.geoPath(mercator, context); const tissot = geoTissot().projection(mercator); const scale = mercator.scale(); const sf = w / 620; const colorScale = d3.scaleSequentialLog(d3.interpolateSpectral).domain([100, 1]); let lonlat = []; for (let lat of d3.range(-80, 81, 10)) { for (let lon of d3.range(-180, 171, 10)) { lonlat.push([lon, lat]); } } function renderMap() { context.clearRect(0, 0, size, size); context.fillStyle = "#dcdcdc"; context.strokeStyle = "#acacac"; context.beginPath(), path({ type: "Sphere" }), context.clip(); context.beginPath(), path(land), context.fill(); context.beginPath(), path(graticule), context.stroke(); tissot({ type: "MultiPoint", coordinates: lonlat}).forEach(indicatrix => { const { coordinates, a, b, theta, s } = indicatrix; const [x, y] = mercator(coordinates); const color = chroma(colorScale(s)); context.save(); context.translate(x, y); context.scale(4 * sf, 4 * sf); context.lineWidth = 1 / 4 * sf; context.fillStyle = color.alpha(0.6); context.strokeStyle = color.alpha(1) context.beginPath(), context.ellipse(0, 0, a, b, theta, 0, 2 * Math.PI); context.fill(); context.stroke(); context.restore(); }); context.fillStyle = "#000"; context.fillText("Aereal scale distortion", size - 105, 10); context.beginPath(), path(sphere), context.stroke(); } for (let i = 0; i < size / 2; ++i) { let s = Math.cosh(i / scale) ** 2, y1 = (i + size / 2 + .5), y2 = size - y1; context.fillStyle = chroma(colorScale(s)).alpha(0.8); context.fillRect(size, y1, 4, 1); context.fillRect(size, y2, 4, 1); } context.fillStyle = "#000"; context.strokeStyle = "#fff"; context.textBaseline = "middle"; [1, 2, 4, 10, 100].forEach((tick, i) => { let y1 = Math.acosh(Math.sqrt(tick)) * scale + size / 2, y2 = size - y1; context.beginPath(); context.moveTo(size, y1); context.lineTo(size + 4, y1); context.stroke(); context.beginPath(); context.moveTo(size, y2); context.lineTo(size + 4, y2); context.stroke(); context.fillText(tick, size + 8, y1 + 1); context.fillText(tick, size + 8, y2 + 1) }) return html`

${d3.select(context.canvas) .call(drag(mercator).on("drag.render", renderMap)) .call(renderMap) .node()}

Areal scale distortion on a Mercator projection

`; }34:T44a,md`
For other projections it might be more useful to calculate *scale distortion* as the ratio of the scale along the meridian or along the parallel at a given point to the scale at a standard point or along a standard line, which is made true to scale. These ratios are called *scale factors*. That along the meridian is called ${tex`h`} and along that the parallel ${tex`k`}. The term *scale error* is frequently applied to ${tex`h-1`} and ${tex`k-1`}. The scale factors ${tex`h`} and ${tex`k`} can also give us some more information about a projection: ${tex`h = k`} if the map is conformal (${tex`\omega = 0`}); ${tex`s = hk`} if meridians and parallels intersect at right angles; ${tex`h = 1 / k`} on an equal-area map if meridians and parallels intersect at right angles.
Let's do a little experiment: What can we find out about the projection of the map below with *Tissot's indicatrix* - without looking at the actual projection definition? `35:T48e,example = { const context = DOM.context2d(w, h); const map = DOM.context2d(w, h); const path = d3.geoPath(projection, map); map.fillStyle = "#dcdcdc"; map.strokeStyle = "#fff"; map.beginPath(), path(land), map.fill(); map.beginPath(), path(countries), map.stroke(); map.strokeStyle = "#acacac"; map.strokeRect(0, 0, w, h); context.fillStyle = "#555"; context.textBaseline = "top"; context.font = "8px monospace"; const format = d3.format("\u2007> 7.2f") const labelWidth = context.measureText("Lon:").width; const lineHeight = 10; const padding = 8; function render(position) { context.clearRect(0, 0, w, h); context.drawImage(map.canvas, 0, 0, w, h); context.fillText("Lon:", padding, padding); context.fillText("Lat:", padding, padding + lineHeight); context.fillText(format(position[0]) + "° W", padding + labelWidth, padding); context.fillText(format(position[1]) + "° N", padding + labelWidth, padding + lineHeight); } render([0, 0]); context.canvas.addEventListener('mousemove', (e) => { render(projection.invert([e.layerX, e.layerY])); }); return context.canvas; }36:T516,md` The [Albers equal-area conic projection](https://en.wikipedia.org/wiki/Albers_projection), here in a specific U.S.-centric configuration. John P. Snyder wrote about it: > *Like other normal conics, the Albers Equal-Area Conic projection has concentric arcs of circles for parallels and equally spaced radii as meridians. [...] The scale factor along meridians is the reciprocal of the scale factor along parallels, to maintain equal area. The standard parallels chosen by Adams for Albers maps of the conterminous United States are lats. 29.5° and 45.5°N. These parallels provide "for a scale error slightly less than 1 percent in the center of the map, with a maximum of 1¼ percent along the northern and southern borders".* That's exactly the result we've got from our little analysis! Not too bad, if I say so to myself. Snyder continues: > *Without measuring the spacing of parallels along a meridian, it is almost impossible to distinguish an unlabeled Albers map of the United States from other conic forms.* But as we have shown, that is possible with the help of *Tissot's indicatrix*. And this is where its real value lies: first of all it is a tool to *analyze* areal and shape distortion, that can be a great help when deciding which projection to choose for a mapping project. `37:Td17,// Formulas are taken from Snyder, Map projections - A working manual, pp. 20-26. // Based on Jason Davies's implementation: // - https://www.jasondavies.com/maps/tissot/ // - https://gis.stackexchange.com/questions/5068/how-to-create-an-accurate-tissot-indicatrix/ geoTissot = { // transpose 2x2 matrix function transpose(mat) { return [mat[0], mat[2], mat[1], mat[3]] } // dot product of two 2x2 matrices function dot(a, b) { return [ a[0] * b[0] + a[1] * b[2], a[0] * b[1] + a[1] * b[3], a[2] * b[0] + a[3] * b[2], a[2] * b[1] + a[3] * b[3] ] } const delta = Math.pow(10, -5.2), rad = Math.PI / 180, fn = function() {}; return function() { let projection = null; function tissot(geojson) { let xy, stream = projection.stream({ point: function(lambda, phi) { xy = [lambda, phi] }, lineStart: fn, lineEnd: fn, polygonStart: fn, polygonEnd: fn }), sf = projection.scale() * rad, // scale factor indicatrices = []; return d3.geoStream(geojson, { point: function(lambda, phi) { const cosPhi = Math.cos(phi * rad); // skip poles (division by zero) if (1e-6 < cosPhi) { let p = (xy = null, stream.point(lambda, phi), xy); if (!!xy) { let dLambda = lambda > 0 ? lambda - delta : lambda + delta, dL = (xy = null, stream.point(dLambda, phi), xy); if (!!xy) { let dPhi = phi > 0 ? phi - delta : phi + delta, dP = (xy = null, stream.point(lambda, dPhi), xy); if (!!xy) { let M = (lambda - dLambda) * sf, m = (phi - dPhi) * sf, // derivates dxdL = (p[0] - dL[0]) / M, dydL = (p[1] - dL[1]) / M, dxdP = (p[0] - dP[0]) / m, dydP = (p[1] - dP[1]) / m, // indicatrix h = Math.sqrt(dxdP * dxdP + dydP * dydP), k = Math.sqrt(dxdL * dxdL + dydL * dydL) / cosPhi, s = -(dydP * dxdL - dxdP * dydL) / cosPhi, A = Math.sqrt(Math.max(0, h * h + k * k + s * 2)), B = Math.sqrt(Math.max(0, h * h + k * k - s * 2)), a = Math.abs((A + B) / 2), b = Math.abs((A - B) / 2), w = Math.asin(Math.max(-1, Math.min(1, B / A))) * 2, // indicatrix rotation mat = dot([dxdL, dxdP, dydL, dydP], [1 / cosPhi, 0, 0, 1]), z = dot(mat, transpose(mat)), theta = Math.atan2(z[2] + z[1], z[0] - z[3]) / 2; indicatrices.push({ coordinates: [lambda, phi], h, k, s, w, a, b, theta }); } } } } } }), indicatrices; } tissot.projection = function(p) { return arguments.length ? (projection = p, tissot) : projection; } return tissot; } }38:T484,contourMapW = { const context = DOM.context2d(w, h); context.textBaseline = "top"; context.font = "10px monospace"; const format = d3.format(".4f"); const labelWidth = context.measureText(`\u03C9 = ${format(0)}`).width; const padding = 8; const map = contourMap(contoursW, magma) context.renderMap = function (indicatrix) { context.clearRect(0, 0, w, h); context.drawImage(map.canvas, 0, 0, w, h); context.fillStyle = "rgba(255, 255, 255, 50%)"; context.fillRect(padding / 2, padding / 2, labelWidth + padding, 16); context.fillStyle = "#333"; context.fillText(`\u03C9 = ${format(indicatrix.w)}`, padding, padding); let [x, y] = projection(indicatrix.coordinates); context.beginPath(), context.arc(x, y, 3, 0, Math.PI * 2), context.stroke(); } let mousedown = false; context.canvas.onmousemove = (e) => { mousedown && pick(e); }; context.canvas.addEventListener("mousedown", (e) => { mousedown = true; context.canvas.onmousemove(e); }); context.canvas.addEventListener("mouseup", (e) => { mousedown = false }); return context; }39:T4a4,contourMapH = { const context = DOM.context2d(w, h); const map = contourMap(contoursH, BrBG); context.textBaseline = "top"; context.font = "10px monospace"; const format = d3.format("7.4f"); const labelWidth = context.measureText(`Scale error (h-1): ${format(0)}`).width; const padding = 8; context.renderMap = function(indicatrix) { context.clearRect(0, 0, w, h); context.drawImage(map.canvas, 0, 0, w, h); context.fillStyle = "rgba(255, 255, 255, 50%)"; context.fillRect(padding / 2, padding / 2, labelWidth + padding, 16); context.fillStyle = "#333"; context.fillText(`Scale error (h-1): ${format(indicatrix.h - 1)}`, padding, padding); let [x, y] = projection(indicatrix.coordinates); context.beginPath(), context.arc(x, y, 3, 0, Math.PI * 2), context.stroke(); } let mousedown = false; context.canvas.onmousemove = (e) => { mousedown && pick(e); }; context.canvas.addEventListener("mousedown", (e) => { mousedown = true; context.canvas.onmousemove(e); }); context.canvas.addEventListener("mouseup", (e) => { mousedown = false; }); return context; }3a:T4a5,contourMapK = { const context = DOM.context2d(w, h); const map = contourMap(contoursK, BrBG); context.textBaseline = "top"; context.font = "10px monospace"; const format = d3.format("7.4f"); const labelWidth = context.measureText(`Scale error (k-1): ${format(0)}`).width; const padding = 8; context.renderMap = function (indicatrix) { context.clearRect(0, 0, w, h); context.drawImage(map.canvas, 0, 0, w, h); context.fillStyle = "rgba(255, 255, 255, 50%)"; context.fillRect(padding / 2, padding / 2, labelWidth + padding, 16); context.fillStyle = "#333"; context.fillText(`Scale error (k-1): ${format(indicatrix.k - 1)}`, padding, padding); let [x, y] = projection(indicatrix.coordinates); context.beginPath(), context.arc(x, y, 3, 0, Math.PI * 2), context.stroke(); } let mousedown = false; context.canvas.onmousemove = (e) => { mousedown && pick(e); }; context.canvas.addEventListener("mousedown", (e) => { mousedown = true; context.canvas.onmousemove(e); }); context.canvas.addEventListener("mouseup", (e) => { mousedown = false; }); return context; }3b:Td4c,function contourMap(contours, colorscale) { const map = DOM.context2d(w, h); const pathA = d3.geoPath(projection, map); const pathB = d3.geoPath(d3.geoIdentity().scale(resolution), map); map.save(); map.beginPath(), map.rect(1, 1, w - 1, h -1 ), map.clip(); map.lineWidth = .5; map.strokeStyle = "#555"; map.setLineDash([5, 2]); contours .forEach((f, i) => { map.save(); if (i > 0) map.globalCompositeOperation = 'source-atop'; map.fillStyle = colorscale(f.value).alpha(0.4); map.beginPath(), pathB(f), map.fill(); map.restore(); //map.beginPath(), pathB(f), map.stroke(); }); map.restore(); let scale = w / (w / resolution + 1); let _w = w / resolution + 1; let _h = h / resolution + 1; for (const contour of contours) { const threshold = contour.value, labels = []; contour.coordinates.forEach(polygon => polygon.forEach((ring, j) => { const p = ring.slice(1, Infinity), // best number of steps to divide ring.length possibilities = d3.range(45, 55), scores = possibilities.map(d => -((p.length - 1) % d)), n = possibilities[d3.scan(scores)], // best starting point: bottom for first rings, top for holes start = 1 + (d3.scan(p.map(xy => (j === 0 ? -1 : 1) * xy[1])) % n), margin = 2; p.forEach((xy, i) => { if ( i % n === start && xy[0] > margin && xy[0] < w - margin && xy[1] > margin && xy[1] < h - margin ) { const a = (i - 2 + p.length) % p.length, b = (i + 2) % p.length, dx = p[b][0] - p[a][0], dy = p[b][1] - p[a][1]; if (dx === 0 && dy === 0) return; labels.push({ threshold, // value xy: xy.map(d => scale * d), angle: Math.atan2(dy, dx), text: `${threshold.toFixed(2)}` }); } }); }) ); // create the mask for this threshold: // the full rectangle minus a rectangle around each label map.save(); map.beginPath(); map.moveTo(0, 0), map.lineTo(w, 0), map.lineTo(w, h), map.lineTo(0, h), map.lineTo(0, 0); const arc = d3.arc(); for (const label of labels) { for (let i = 0; i < 2 * Math.PI; i += 0.2) { const pos = [Math.cos(i) * 13, -Math.sin(i) * 10], c = Math.cos(label.angle), s = Math.sin(label.angle); map[i === 0 ? "moveTo" : "lineTo"]( label.xy[0] + pos[0] * c - pos[1] * s, label.xy[1] + pos[1] * c + pos[0] * s ); } } //map.stroke(); // uncomment to see the mask map.clip(); map.lineWidth = .5; map.strokeStyle = "#555"; map.setLineDash([5, 2]); map.beginPath(), pathB(contour), map.stroke(); // draw labels for this threshold map.restore(); for (const label of labels) { addlabel(map, label); } } map.globalCompositeOperation = 'multiply'; map.strokeStyle = "#fff"; map.beginPath(), pathA(countries), map.stroke(); map.globalCompositeOperation = 'destination-over'; map.fillStyle = "#dcdcdc"; map.beginPath(), pathA(land), map.fill(); map.globalCompositeOperation = 'source-over'; map.strokeStyle = "#acacac"; map.strokeRect(0, 0, w, h); return map; }3c:[] 2e:{"title":"Tissot's indicatrix","cells":[{"id":3911,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":3540,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":3892,"value":"viewof lon = slider({\n min: -180, \n max: 180, \n step: 0.1,\n value: 30,\n description: \"Longitude\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":3894,"value":"viewof lat = slider({\n min: -90, \n max: 90, \n step: 0.1,\n value: 45,\n description: \"Latitude\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":605,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1279,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":1767,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":1901,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":1935,"value":"md`To start our investigation we create a set of control points, that we can pass to the \\`geoTissot\\` generator function. It's a regular grid of points over the canvas area, reprojected to longitude and latitude and formatted as a GeoJSON object.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1981,"value":"geojson = {\n let grid = [];\n for (let y of d3.range(0, h + resolution, resolution)) {\n for (let x of d3.range(0, w + resolution, resolution)) {\n grid.push(projection.invert([x, y]));\n }\n }\n return { type: \"MultiPoint\", coordinates: grid };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2342,"value":"md`The next step is to initialize the generator function and specify the [D3 projection](https://github.com/d3/d3-geo#projections) for which we want to create the *indicatrices*:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1991,"value":"tissot = geoTissot().projection(projection)","pinned":true,"hidden":false,"mode":"ojs"},{"id":3374,"value":"md`Finally we pass the grid of control points to the indicatrix generator. It accepts any GeoJSON object and computes indicatrices for every coordinate pair of the input geometry.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2352,"value":"data = tissot(geojson);","pinned":true,"hidden":false,"mode":"ojs"},{"id":3386,"value":"md`The return value is an array of *indicatrices* with the following properties:\n\n\\`coordinates\\` - longitude and latitude of the indicatrix \n\\`h\\` - scale factor along the meridian ${tex`h`} \n\\`k\\` - scale factor along the parallel ${tex`k`} \n\\`w\\` - maximum angular distortion ${tex`\\omega`} \n\\`s\\` - areal scale distortion ${tex`s`} \n\\`a\\` - length of semi-major axis ${tex`a`} \n\\`b\\` - length of semi-minor axis ${tex`b`} \n\\`theta\\` - indicatrix rotation in radians\n\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1948,"value":"{\n const context = DOM.context2d(w, h);\n const path = d3.geoPath(projection, context).pointRadius(3);\n \n context.fillStyle = \"#dcdcdc\";\n context.strokeStyle = \"#fff\";\n \n context.beginPath(), path(land), context.fill();\n context.beginPath(), path(countries), context.stroke();\n context.strokeStyle = \"#acacac\";\n context.strokeRect(0, 0, w, h);\n \n context.beginPath();\n context.rect(1, 1, w - 2, h - 2);\n context.clip()\n \n context.fillStyle = \"#555\";\n data.forEach(indicatrix => {\n const { coordinates, a, b, theta } = indicatrix;\n const [x, y] = projection(coordinates);\n context.beginPath(), context.ellipse(x, y, a, b, theta, 0, Math.PI * 2), context.fill();\n })\n \n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2354,"value":"md`The mapped indicatrices don't look that spectacular at first sight. In fact, it wouldn't make much sense to show magnified indicatrices on a map at this scale. The scale errors are so small that differences would barely be visible and the ellipses would more or less appear as circles. However, we can now start analyzing the generated data.\n\nOn a *[conformal](https://en.wikipedia.org/wiki/Conformal_map)* map projection there is no local shape distortion. Thus, any small figure on the earth is nearly similar to its image on the map, but areas are not preserved. Altough this is a very nice property for a projection, this projection doesn't seem to be conformal:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2356,"value":"conformal = {\n return data\n .map(indicatrix => indicatrix.w)\n .every(omega => omega < 1e-6)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2385,"value":"md`An *equal-area* map projection, on the other hand, shows areas on the map proportional to their area on the globe. It generally distorts shapes in order to do that, so no projection can be both equal-area and conformal. A projection is of course not necessarily equal-area only because it isn't conformal, but in this case it is true:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2387,"value":"equalarea = {\n return data\n .map(indicatrix => Math.abs(indicatrix.s - 1))\n .every(delta => delta < 1e-6)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2402,"value":"md`Altough equal-area projections distort angles on a map, that doesn't mean that *all* angles except the angle bewteen the semi-major and semi-minor axis of the ellipse are distorted. On a conformal projection, parallels and meridians always cross rectangularly on the map. And this is also the case here:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2404,"value":"rightAngles = {\n return data\n .map(indicatrix => Math.abs(indicatrix.s - indicatrix.h * indicatrix.k))\n .every(delta => delta < 1e-6)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":3241,"value":"md`
Now that we know we have an equal-area projection with meridians and parallels intersecting at right angles, let's examine the distortion of shapes in more detail. A nice way of doing it is to visualize the *maximum angular distortion* ${tex`{\\omega}`} with [contour lines](https://en.wikipedia.org/wiki/Contour_line) on the map.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2055,"value":"contoursW = d3.contours()\n .thresholds([0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06])\n .size([w / resolution + 1, h / resolution + 1])\n (data.map(indicatrix => indicatrix.w));","pinned":false,"hidden":false,"mode":"ojs"},{"id":4303,"value":"contourMapW.canvas","pinned":false,"hidden":false,"mode":"ojs"},{"id":2809,"value":"md`Visible are two bands streching from west to east, where the angular distortion stays below 0.01. While the maximum angular distortion quickly increases in the northern and southern regions of the map, it stays below a threshold value of 0.02 in the central area.\n\nAt this point of time it's probably safe to say that it is a [conic projection](https://en.wikipedia.org/wiki/Map_projection#Conic) using two standard parallels, where areas *and* angles are preserved. Nonetheless one might be further interested in *how* angles are distorted on the map and examine the *scale errors* along the meridians and the parallels.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2812,"value":"contoursH = d3.contours()\n .thresholds([-0.04, -0.03, -0.02, -0.01, 0, 0.01])\n .size([w / resolution + 1, h / resolution + 1])\n (data.map(indicatrix => indicatrix.h - 1));","pinned":false,"hidden":false,"mode":"ojs"},{"id":3130,"value":"contoursK = d3.contours()\n .thresholds([-0.01, 0, 0.01, 0.02, 0.03, 0.04])\n .size([w / resolution + 1, h / resolution + 1])\n (data.map(indicatrix => indicatrix.k - 1));","pinned":false,"hidden":false,"mode":"ojs"},{"id":4318,"value":"contourMapH.canvas","pinned":false,"hidden":false,"mode":"ojs"},{"id":4326,"value":"contourMapK.canvas","pinned":false,"hidden":false,"mode":"ojs"},{"id":3285,"value":"md`The lines of true scale are now visible in both maps (threshold value 0) and seperate the northern and southern regions from the central area of the map, where the scale error is below 1% for both, along the meridians and the parallels. In the northern and southern regions of the map shapes are stretched along the parallel and contracted along the meridian. The opposite is true for the central area of the map, but to a much lesser extend.\n\nAs a final step in our investigation let's make a guess where the standard parallels of the projection are located approximately:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3187,"value":"parallels = new Set(data\n .filter(indicatrix => indicatrix.w < 1e-5)\n .map(indicatrix => Math.round(indicatrix.coordinates[1] * 1000) / 1000)\n );","pinned":true,"hidden":false,"mode":"ojs"},{"id":3309,"value":"md`The projection is thus a conic equal area projection that preserves shape and area along the standard parallels at approximately 29.5° and 45.5° N. Meridians and parallels intersect at right angles. The scale error in the central area of the map is kept below 1% along the meridians and the parallels. And the projection is:\n\n(drumroll) \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1899,"value":"projection = d3.geoAlbers()\n .scale(640 * w / 640)\n .translate([w / 2, h / 2]);","pinned":true,"hidden":false,"mode":"ojs"},{"id":4382,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":1879,"value":"md`

References: \n\n\n*John P. Snyder, *[Map projections - A working manual](https://pubs.er.usgs.gov/publication/pp1395)*, Washington 1987, pp. 20-26, 98-103. \nJason Davies:* [Tissot’s Indicatrix](https://www.jasondavies.com/maps/tissot/)*, retrieved: 24/05/2019. \nStackoverflow:* [How to create an accurate Tissot Indicatrix?](https://gis.stackexchange.com/questions/5068/how-to-create-an-accurate-tissot-indicatrix/)*, retrieved: 22/05/2019.*\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3280,"value":"md`---\n## Implementation`","pinned":false,"hidden":false,"mode":"ojs"},{"id":943,"value":"$37","pinned":true,"hidden":false,"mode":"ojs"},{"id":1460,"value":"md`---\n## Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1904,"value":"w = Math.min(640, width);","pinned":false,"hidden":false,"mode":"ojs"},{"id":1907,"value":"h = Math.floor(w / 4 * 2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2659,"value":"resolution = 5","pinned":false,"hidden":false,"mode":"ojs"},{"id":4188,"value":"magma = chroma.scale([\"#fcfdbf\", \"#b73779\", \"#000004\"]).domain([0, 0.05]).mode('lch');","pinned":false,"hidden":false,"mode":"ojs"},{"id":4196,"value":"BrBG = chroma.scale([\"#543005\", \"#f5f5f5\", \"#003c30\"]).domain([-0.04, 0, 0.03]);","pinned":false,"hidden":false,"mode":"ojs"},{"id":4422,"value":"function addlabel(context, label) {\n context.save();\n context.translate(...label.xy);\n context.rotate(label.angle + (Math.cos(label.angle) < 0 ? Math.PI : 0));\n context.textAlign = \"center\";\n context.fillStyle = \"#555\";\n context.font = \"8px sans-serif\";\n context.fillText(label.text, 0, 5);\n context.restore();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2080,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":2816,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":3132,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":4169,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":4296,"value":"update = {\n let indicatrix = tissot({ type: \"Point\", coordinates: position })[0];\n contourMapW.renderMap(indicatrix);\n contourMapH.renderMap(indicatrix);\n contourMapK.renderMap(indicatrix);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4308,"value":"viewof position = new View([-96.5, 38.5])","pinned":false,"hidden":false,"mode":"ojs"},{"id":3441,"value":"function pick(event) {\n viewof position.value = projection.invert([event.layerX, event.layerY]);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1305,"value":"function drag(projection) {\n let v0, q0, r0;\n \n function dragstarted() {\n v0 = versor.cartesian(projection.invert(d3.mouse(this)));\n q0 = versor(r0 = projection.rotate());\n }\n \n function dragged() {\n const v1 = versor.cartesian(projection.rotate(r0).invert(d3.mouse(this)));\n const q1 = versor.multiply(q0, versor.delta(v0, v1));\n projection.rotate(versor.rotation(q1));\n }\n \n return d3.drag()\n .on(\"start\", dragstarted)\n .on(\"drag\", dragged);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3805,"value":"rad = Math.PI / 180","pinned":false,"hidden":false,"mode":"ojs"},{"id":3808,"value":"deg = 180 / Math.PI","pinned":false,"hidden":false,"mode":"ojs"},{"id":1372,"value":"sphere = ({ type: \"Sphere\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"graticule = d3.geoGraticule10()","pinned":false,"hidden":false,"mode":"ojs"},{"id":1920,"value":"countries = topojson.mesh(world, world.objects.countries, (a, b) => a !== b)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1917,"value":"land = topojson.feature(world, world.objects.land)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3302,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"world = fetch(\"https://unpkg.com/world-atlas@1/world/110m.json\")\n .then(response => response.json())","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"d3 = require(\"d3@5\", \"d3-geo-projection\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"topojson = require(\"topojson-client@3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1307,"value":"versor = require(\"versor@0.0.3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1525,"value":"chroma = require(\"chroma-js\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":3829,"value":"import {slider} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":4305,"value":"import {View} from \"@mbostock/synchronized-views\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3c","stdlib":"1","titling":"auto"}