CartoLab Academic Reference Manual

Mathematical formulations, theoretical foundations, and verified literature for all 13 Processing algorithms and 5 advanced Layout Studio systems. Version 2.8.6.

This manual provides the complete theoretical apparatus behind every algorithm and layout engine in PlanX CartoLab. Each entry supplies the mathematical derivation, parameter semantics, interpretative guidance, and peer-reviewed references. All equations are typeset with MathJax 3.

geometric_interval_classification Classification

1. Advanced Classification (GIC / Head-Tail / Fisher-Jenks)

Theoretical Background

Classification is the fundamental cartographic operation that assigns continuous data values to discrete visual classes. The choice of classification method profoundly influences the map reader's perception of spatial pattern. Traditional equal-interval or quantile methods often fail for urban datasets that exhibit strong positive skew, power-law tails, or multimodal distributions. This algorithm implements three complementary approaches that together cover the principal distributional forms encountered in planning and demographic data.

The Adaptive Geometric Interval Classification (GIC) fits a geometric series $C_k = a \cdot r^k$ to the data, optimising the ratio $r$ and scale $a$ to minimise the sum of squared deviations between observed class frequencies and the expected mean frequency $\bar{f} = n/K$. The method performs a coarse-then-fine grid search over feasible $r$ values, using a shift parameter to handle datasets containing zero or negative values. GIC is particularly effective for variables such as population density, land value, and building height, where the distribution is approximately log-normal.

Head/Tail Breaks (Jiang, 2013) recursively partitions the data at the arithmetic mean, retaining the head (values above the mean) for further subdivision. The recursion terminates when the head is no longer "heavy-tailed"—defined as fewer than 2 items or less than 40% (head proportion ≥ 80%) of the current subset. This method is designed for data with a heavy-tailed or power-law distribution, common in city size hierarchies, traffic flows, and income data. Unlike arbitrary class counts, Head/Tail Breaks determines the number of classes from the data structure. Fisher-Jenks Natural Breaks (Jenks, 1967; Fisher, 1958) uses dynamic programming to minimise the sum of within-class variances. For $n > 5000$ observations, the algorithm automatically samples to maintain performance while preserving distributional fidelity.

Mathematical Formulation

Geometric Interval (GIC) objective. $$ \mathcal{L}(r, a) = \sum_{j=1}^{K} \left( f_j - \frac{n}{K} \right)^2, \qquad C_k = a \cdot r^k $$ where $f_j$ is the number of observations falling in class $j$, $n$ is the total count, and $K$ is the number of classes. The pair $(r^*, a^*)$ that minimises $\mathcal{L}$ is found by coarse search over $r \in [0.5, 2.5]$ followed by fine-grained refinement. \tag{1}
Head/Tail Breaks mean partition. $$ \mu = \frac{1}{n} \sum_{i=1}^{n} x_i, \qquad H = \{\,x_i \mid x_i > \mu\,\}, \qquad T = \{\,x_i \mid x_i \leq \mu\,\} $$ Recursion continues on $H$ while $|H|/|H \cup T| \geq 0.8$ and $|H| \geq 2$. \tag{2}
Fisher-Jenks dynamic programming recurrence. $$ \mathrm{SSQ}(i, j) = \sum_{k=i}^{j-1} (x_k - \bar{x}_{ij})^2, \qquad \bar{x}_{ij} = \frac{1}{j-i} \sum_{k=i}^{j-1} x_k $$ $$ \mathrm{DP}[k][j] = \min_{i \,<\, j} \bigl( \mathrm{DP}[k-1][i] + \mathrm{SSQ}(i, j) \bigr) $$ where $\mathrm{DP}[k][j]$ is the minimum within-class variance when partitioning the first $j$ sorted observations into $k$ classes. Computed via cumulative-sum acceleration for $\mathcal{O}(K n^2)$ time. \tag{3}
Within-class variance minimisation (Jenks optimality criterion). $$ \mathrm{GVF} = 1 - \frac{\sum_{c=1}^{K} \sum_{x_i \in C_c} (x_i - \bar{x}_c)^2} {\sum_{i=1}^{n} (x_i - \bar{x})^2} $$ where $\mathrm{GVF}$ is the Goodness of Variance Fit (0–1), higher is better. Jenks breaks maximise GVF for a given $K$. \tag{4}
Geometric interval with non-positive shift. $$ x_i' = x_i + |\min(0, x_{\min})| + \varepsilon, \qquad \varepsilon = 10^{-6} $$ Shift is applied transparently; class boundaries are mapped back to the original scale before output. \tag{5}
Maximum Distribution Breaks (Largest Gap Partition). $$ \Delta_i = x_{(i+1)} - x_{(i)}, \qquad \text{for sorted } x_{(1)} \le x_{(2)} \le \dots \le x_{(n)} $$ Partitions the data at the $K-1$ largest spacing differentials $\Delta_i$, isolating natural clusters. \tag{6}
Pretty Nice-Round Breaks (Heckbert Algorithm). $$ \Delta = \frac{x_{\max} - x_{\min}}{K}, \qquad d = 10^{\lfloor \log_{10} \Delta \rfloor} $$ $$ \text{Step } s \in \{1, 2, 5, 10\} \times d \quad \text{such that } s \approx \Delta $$ Yields aesthetically intuitive, human-readable integer class limits ($10, 20, 50, 100, 250\dots$). \tag{7}

Parameters

ParameterTypeDefaultDescription
Input layerVector (any geometry)—Layer containing the field to classify
Field to classifyNumeric field—Must be integer or real; NULL values are excluded
Number of classesInteger5Range 2–20; Head/Tail Breaks ignores this and auto-determines
Classification methodEnumGIC0 = Geometric Interval, 1 = Head/Tail Breaks, 2 = Fisher-Jenks, 3 = Maximum Breaks, 4 = Pretty Breaks
Color rampEnumViridisViridis, Plasma, Inferno, Magma, Cividis, Turbo, Mako, Rocket, Blues, Oranges, YlOrRd, Purples, Greens

Output Description

The output layer is a copy of the input with two appended fields: gic_class (integer, 0-indexed class assignment) and gic_label (string, formatted range label such as “123.40 – 567.80”). The layer is automatically styled with a graduated renderer using the chosen colour ramp and clean class-boundary legend entries.

Interpretation Guide

Interpretation Choose GIC for continuous skewed data (population, density, value), Head/Tail Breaks when the histogram shows a heavy right tail (city sizes, traffic counts, income), Fisher-Jenks for general-purpose choropleth mapping where visual homogeneity within classes is the priority, Maximum Breaks when highlighting distinct multimodal clusters, and Pretty Breaks for public-facing maps requiring intuitive, rounded interval boundaries. Examine the GVF value reported in the log: values above 0.8 indicate a strong class structure.

References

  1. ESRI. Geometric Interval Classification. ArcGIS Desktop Help: Classification methods. ESRI Technical Paper, 2005.
  2. Jenks, G. F. The Data Model Concept in Statistical Mapping. International Yearbook of Cartography, 7: 186–190, 1967.
  3. Fisher, W. D. On Grouping for Maximum Homogeneity. Journal of the American Statistical Association, 53(284): 789–798, 1958. doi:10.1080/01621459.1958.10501479
  4. Jiang, B. Head/Tail Breaks: A New Classification Scheme for Data with a Heavy-Tailed Distribution. The Professional Geographer, 65(3): 482–494, 2013. doi:10.1080/00330124.2012.700499
  5. Heckbert, P. S. Nice Numbers for Graph Labels. Graphics Gems, Academic Press, pp. 61–63, 1990.
  6. Slocum, T. A., McMaster, R. B., Kessler, F. C., & Howard, H. H. Thematic Cartography and Geovisualization (3rd ed.). Pearson, 2009. ISBN 978-0132298346.
bivariate_choropleth Thematic Mapping

2. Bivariate Choropleth Map

Theoretical Background

Bivariate choropleth maps encode the relationship between two spatially coincident variables within a single areal symbolisation. Unlike side-by-side univariate maps, a bivariate matrix forces the reader to perceive co-variation as a single visual gestalt. The technique was formalised by the U.S. Census Bureau in the 1970s and refined by Olson (1981), who demonstrated that spectrally encoded two-variable maps could be read accurately when the colour space is perceptually orthogonal.

The CartoLab implementation classifies each of two numeric fields independently into $N$ bins using the selected method (Geometric Interval or Fisher-Jenks), producing an $N \times N$ grid of cross-classified observations. The colour for each cell $(r, c)$ is computed via bilinear interpolation in RGB space between four user-specified corner colours: Low-Low $(C_{LL})$, Low-High $(C_{LH})$, High-Low $(C_{HL})$, and High-High $(C_{HH})$. This interpolation produces a perceptually smooth colour surface that preserves the logical ordering of both variables. The output layer receives an auto-generated categorised renderer with $N^2$ legend entries, each displaying a swatch from the colour matrix.

The theoretical underpinning draws on colour-order systems developed by Eyton (1984) for complementary-colour maps and Leonowicz (2006) for two-variable choropleth cartography. Care is taken to ensure the diagonal (equal-class cells) is visually coherent and that off-diagonal cells maintain sufficient discriminability. For $N \geq 5$, the matrix produces 25 classes, at which point individual cell discrimination becomes challenging; the recommended grid size is 3 or 4.

Mathematical Formulation

Bilinear colour interpolation in RGB space. $$ \mathbf{C}(r, c) = (1 - \hat{y})(1 - \hat{x})\,\mathbf{C}_{LL} + (1 - \hat{y})\,\hat{x}\,\mathbf{C}_{LH} + \hat{y}(1 - \hat{x})\,\mathbf{C}_{HL} + \hat{y}\,\hat{x}\,\mathbf{C}_{HH} $$ where $\hat{x} = c/(N-1)$ and $\hat{y} = r/(N-1)$ are the normalised class indices (0 to 1), and each $\mathbf{C}$ is an RGB vector. \tag{1}
Independently classified axes. $$ x_i^{\text{class}} = \operatorname{classify}(x_i; N, \text{method}), \qquad y_i^{\text{class}} = \operatorname{classify}(y_i; N, \text{method}) $$ The composite bivariate class index is $b_i = x_i^{\text{class}} \cdot N + y_i^{\text{class}}$ (0-indexed). \tag{2}
Class boundary recovery. For each axis, class breaks $\{b_0, b_1, \ldots, b_N\}$ are computed via the selected classification method. A feature is assigned to class $k$ if $b_k \leq x_i < b_{k+1}$. \tag{3}
Matrix cell assignment. $$ M[r][c] = \{\,i \mid x_i \in \text{class } c,\; y_i \in \text{class } r\,\} $$ where $r, c \in \{0, \ldots, N-1\}$ and $M[r][c]$ is the set of feature indices in that matrix cell. \tag{4}

Parameters

ParameterTypeDefaultDescription
Input layerPolygon vector—Contains both variables to cross-classify
X-axis variableNumeric field—Horizontal axis in the colour matrix
Y-axis variableNumeric field—Vertical axis in the colour matrix
Grid sizeInteger3$N$ (2–7); produces $N^2$ classes
Classification methodEnumGeometric IntervalApplied independently to both variables
4 corner coloursColour (RGBA)teal/gold/blue/copperUser-adjustable LL, LH, HL, HH corner anchors

Output Description

The output polygon layer carries three new fields: bivar_x_class (integer, 0-indexed X-axis bin), bivar_y_class (integer, 0-indexed Y-axis bin), and bivar_class (integer, 0 to $N^2-1$, the composite colour-matrix cell index). The layer is auto-styled with a categorised renderer using the $N^2$ bilinearly interpolated colours.

Interpretation Guide

Interpretation Read the matrix diagonally (bottom-left to top-right): features on the diagonal have concordant rankings on both variables. Off-diagonal cells indicate discordance. The LL corner (bottom-left) represents low values on both variables; HH (top-right) is high on both. Use $N = 3$ for a clean, readable 9-class legend; $N = 4$ for finer discrimination in a publication figure. Always pair with a bivariate legend (CartoLab dashboard can generate one) because the colour semantics are not self-evident to map readers unfamiliar with the technique.

References

  1. Olson, J. M. Spectrally Encoded Two-Variable Maps. Annals of the Association of American Geographers, 71(2): 259–276, 1981. doi:10.1111/j.1467-8306.1981.tb01352.x
  2. Eyton, J. R. Complementary-Color Two-Variable Maps. Annals of the Association of American Geographers, 74(3): 477–490, 1984.
  3. Leonowicz, A. Two-Variable Choropleth Maps as a Useful Tool for Visualization of Geographical Relationship. Geografija, 42(1): 33–37, 2006.
  4. Brewer, C. A. Designing Better Maps: A Guide for GIS Users (2nd ed.). ESRI Press, 2015. ISBN 978-1589484405.
  5. Carstensen, L. W. Bivariate Choropleth Mapping: The Effects of Axis Scaling. The American Cartographer, 13(1): 27–41, 1986.
  6. Dunn, R. A Dynamic Approach to Two-Variable Color Mapping. The American Statistician, 43(4): 245–252, 1989.
  7. Elmer, M. E. Symbol Considerations for Bivariate Thematic Maps. Cartography and Geographic Information Science, 39(3): 140–150, 2012.
hexbin_aggregate Aggregation

3. Hexbin Aggregation

Theoretical Background

Hexagonal binning is a spatial aggregation technique that partitions a two-dimensional point set into a regular tessellation of regular hexagons. The method was popularised in statistical graphics by Carr et al. (1987) as a solution to overplotting in large scatterplots and later extended to geographical space. Hexagons offer several advantages over square grids: each cell has six equidistant neighbours (versus four orthogonal plus four diagonal for squares), producing smoother density estimates; the circular packing of hexagons reduces sampling bias from edge effects; and the regular tessellation avoids the directional artefacts that rectangular grids can introduce.

The CartoLab implementation uses a pointy-top hexagonal lattice in axial coordinates $(q, r)$, derived from the standard cube-coordinate system $(x, y, z)$ where $x + y + z = 0$. The conversion from geographic coordinates to axial coordinates uses the hexagon circumradius $R$ (distance from centre to vertex). Cube rounding ensures that each point is assigned to the nearest hexagon centre. Unlike some implementations that emit a full grid, CartoLab only emits cells that contain at least one point, conserving output size for sparse datasets.

The method aggregates a user-selectable statistic: count (the number of points falling in each cell), sum (the total of an optional weight field), or mean (sum divided by count). The output hexagons are rendered as flat-topped (rotated 90° from the internal pointy-top computation) for consistency with common cartographic convention. Each hexagon is represented as a regular polygon with six vertices computed at 60° intervals.

Mathematical Formulation

Point to axial coordinates (pointy-top). $$ q = \frac{\frac{\sqrt{3}}{3}\,x - \frac{1}{3}\,y}{R}, \qquad r = \frac{\frac{2}{3}\,y}{R} $$ where $(x, y)$ are map coordinates relative to the layer CRS origin and $R$ is the hexagon circumradius. \tag{1}
Cube rounding to nearest hex cell. $$ x_q = q,\quad x_r = r,\quad x_s = -q - r $$ Round each coordinate: $\hat{x}_q = \lfloor x_q \rceil$, etc. If $\hat{x}_q + \hat{x}_r + \hat{x}_s \neq 0$, correct the coordinate with the largest rounding residual: $$ \text{reassign the component with max } |x_k - \hat{x}_k| \text{ to restore } \hat{x}_q + \hat{x}_r + \hat{x}_s = 0 $$ Final axial: $(\hat{q}, \hat{r}) = (\hat{x}_q, \hat{x}_r)$. \tag{2}
Hexagon centre from axial coordinates (pointy-top). $$ x_{\text{ctr}} = R \cdot \left( \sqrt{3}\,q + \frac{\sqrt{3}}{2}\,r \right), \qquad y_{\text{ctr}} = R \cdot \left( \frac{3}{2}\,r \right) $$ \tag{3}
Vertex positions at 60° intervals (pointy-top). $$ v_k = \left(x_{\text{ctr}} + R\cos\theta_k,\; y_{\text{ctr}} + R\sin\theta_k\right), \qquad \theta_k = \frac{\pi}{6} + \frac{k\pi}{3},\; k = 0, \ldots, 5 $$ \tag{4}
Aggregation statistics per cell $c$. $$ n_c = |\{i \mid \operatorname{cell}(p_i) = c\}|, \qquad s_c = \sum_{\operatorname{cell}(p_i) = c} w_i, \qquad \bar{w}_c = s_c / n_c $$ \tag{5}

Parameters

ParameterTypeDefaultDescription
Input point layerPoint vector—Points to bin
Hexagon radiusFloat1000Circumradius in map units (CRS-dependent)
Weight fieldNumeric field (optional)—Used for sum and mean statistics
StatisticEnumCount0 = Count, 1 = Sum, 2 = Mean

Output Description

The output is a polygon layer of regular hexagons. Non-empty cells carry hex_q and hex_r (axial coordinates), hex_count (integer, always populated), hex_sum (weight-field sum, present if sum or mean selected), and hex_mean (present if mean selected). The layer is auto-graduated with the Viridis colour ramp on the active statistic field.

Interpretation Guide

Interpretation Choose the radius to balance spatial resolution against visual clarity: a good starting point is 1/30th to 1/50th of the map extent diagonal. The hexbin map reveals spatial density gradients that point maps obscure. Hex grids with a Viridis ramp produce continuous-looking surfaces suitable for publication. Because only occupied cells are emitted, large gaps in the grid indicate genuinely empty regions, not missing data.

References

  1. Carr, D. B., Littlefield, R. J., Nicholson, W. L., & Littlefield, J. S. Scatterplot Matrix Techniques for Large N. Journal of the American Statistical Association, 82(398): 424–436, 1987. doi:10.1080/01621459.1987.10478445
  2. Carr, D. B., Olsen, A. R., & White, D. Hexagon Mosaic Maps for Display of Univariate and Bivariate Geographical Data. Cartography and Geographic Information Systems, 19(4): 228–236, 1992. doi:10.1559/152304092783721231
  3. Birch, C. P. D., Oom, S. P., & Beecham, J. A. Rectangular and Hexagonal Grids Used for Observation, Experiment and Simulation in Ecology. Ecological Modelling, 206(3–4): 347–359, 2007.
  4. Battersby, S. E., Strebe, D., & Finn, M. P. Shapes on a Plane: Evaluating Hexagonal Binning and Other Tessellation Strategies. Cartographic Perspectives, 85: 27–40, 2017.
  5. Lewin-Koh, N. Hexagon Binning: An Overview. R package hexbin vignette, 2011.
  6. Red Blob Games (Patel, A.). Hexagonal Grids. 2013 (revised 2023). Interactive exposition of axial/cube coordinate systems.
  7. Scott, D. W. Multivariate Density Estimation: Theory, Practice, and Visualization (2nd ed.). Wiley, 2015. ISBN 978-0471697558.
dot_density Thematic Mapping

4. Dot-Density Map

Theoretical Background

The dot-density map is one of the oldest thematic mapping techniques, dating to the earliest demographic atlases of the 19th century. Each dot represents a fixed number of units (e.g., one dot = 500 persons), and dots are placed within enumeration areas to convey both the magnitude and the approximate spatial distribution of a phenomenon. Unlike a choropleth map, which forces the reader to infer density from colour intensity across administrative boundaries, a dot-density map shows within-unit variation that can reveal settlement patterns, ecological distributions, and other spatially heterogeneous phenomena.

The CartoLab implementation generates dots via seeded random stratified sampling within each polygon's bounding box, rejecting dots that fall outside the polygon (hole-aware with OGC-compliant point-in-polygon testing). The number of dots per feature is $n_{\text{dots}} = \lfloor v_i / d \rfloor$ where $v_i$ is the count field value and $d$ is the user-specified value per dot. Each feature uses a deterministic pseudo-random seed derived from the user seed and the feature's row index ($\text{seed}_i = \text{seed} \times 1{,}000{,}003 + i$), ensuring reproducible dot placement across runs.

Robinson et al. (1984) recommend limiting dot density to approximately 200 dots/cm² in print to avoid coalescence; for screen display, a lower threshold applies. CartoLab scales $d$ automatically if $n_{\text{dots}}$ would exceed a performance ceiling, and reports the adjustment. Lavin (1986) explored the perceptual accuracy of dot-density maps and found that random placement within enumeration areas is visually acceptable when the dot value is small relative to the feature size, which is the typical use case.

Mathematical Formulation

Number of dots per feature. $$ n_i^{\text{dots}} = \left\lfloor \frac{v_i}{d} \right\rfloor $$ where $v_i$ is the count-field value for feature $i$ and $d$ is the user-specified value per dot. Remainder rounding is truncated. \tag{1}
Seeded deterministic random placement. $$ \text{seed}_i = S \times 1{,}000{,}003 + i $$ where $S$ is the user seed and $i$ is the 0-indexed feature row number. This produces a repeatable pseudo-random sequence per feature across runs. \tag{2}
Bounding-box stratified sampling. For each of $n_i^{\text{dots}}$ trials, generate a candidate point $$ p_{\text{cand}} = \bigl( x_{\min} + u \cdot w,\; y_{\min} + v \cdot h \bigr) $$ where $u, v \sim \mathcal{U}(0, 1)$ (derived from the PRNG), $w = x_{\max} - x_{\min}$, $h = y_{\max} - y_{\min}$. Accept $p_{\text{cand}}$ iff $\operatorname{contains}(\text{geom}_i, p_{\text{cand}}) = \text{true}$. \tag{3}
Point-in-polygon (ray-casting). For a candidate point $p$, cast a horizontal ray to $+\infty$. Count intersections with polygon edges. $p$ is inside iff the count is odd, with special handling for edge-on-vertex and colinear edge cases per the even-odd rule (OGC Simple Features specification). \tag{4}

Parameters

ParameterTypeDefaultDescription
Input polygon layerPolygon vector—Enumeration areas containing the count field
Count fieldNumeric field—The magnitude to represent as dots
Value per dotFloat1Each dot represents this many units
Random seedInteger42Deterministic seed for reproducible placement

Output Description

The output is a point layer where each dot is a single point geometry. Each dot inherits all attributes from its source polygon, enabling multi-group dot maps (e.g., race/ethnicity categories) by filtering on a category field after generation. No new fields are added to the dot layer beyond the source attributes.

Interpretation Guide

Interpretation Dot-density maps excel at conveying the texture of a distribution: clusters, voids, and gradients are visible at a glance. Choose the dot value such that the densest polygon receives roughly 20–50 dots, ensuring individual dots remain distinguishable. For multi-group maps, generate separate dot layers and stack them, using distinct hues. Always note the dot value prominently in the legend, as readers have no visual anchor for absolute magnitudes without it.

References

  1. Robinson, A. H., Sale, R. D., Morrison, J. L., & Muehrcke, P. C. Elements of Cartography (6th ed.). Wiley, 1995. ISBN 978-0471555797.
  2. Lavin, S. J. Mapping Continuous Geographical Distributions Using Dot-Density Shading. The American Cartographer, 13(2): 140–150, 1986.
  3. Dent, B. D., Torguson, J. S., & Hodler, T. W. Cartography: Thematic Map Design (6th ed.). McGraw-Hill, 2009. ISBN 978-0072943825.
  4. Mackay, J. R. An Analysis of Isopleth and Dot Maps. Annals of the Association of American Geographers, 45(3): 296–297, 1955.
  5. Heywood, I., Cornelius, S., & Carver, S. An Introduction to Geographical Information Systems (4th ed.). Pearson, 2011. ISBN 978-0273722595.
  6. Krygier, J. & Wood, D. Making Maps: A Visual Guide to Map Design for GIS (3rd ed.). Guilford Press, 2016. ISBN 978-1462509980.
compute_cartogram Cartogram

5. Continuous-Area Cartogram (Diffusion)

Theoretical Background

A continuous-area cartogram distorts the geometry of enumeration units so that each unit's area is proportional to a selected variable (e.g., population, GDP) while preserving, as far as possible, the contiguity and shape of the original map. The technique has a long pre-computational history but was rendered practical for arbitrary datasets by Gastner and Newman's (2004) diffusion method, which treats the map as a physical medium in which a "density" flows from regions of excess area to regions of deficit until the area error is eliminated.

The CartoLab implementation follows the Gastner–Newman algorithm closely. Each polygon $i$ has a target area $A_i^{\text{target}}$ proportional to its variable value. The current "mass" $m_i$ is defined as the discrepancy between the current equivalent radius and the target radius. A Gaussian-weighted force kernel distributes this mass to neighbouring boundaries, displacing vertices in the direction that reduces area discrepancy. The process iterates until either the mean size error (the ratio of max to min of current and target area) falls below a threshold or the maximum iteration count is reached.

The algorithm employs a spatial hash grid for $\mathcal{O}(n)$ neighbour lookup, replacing the naïve $\mathcal{O}(n^2)$ search. Displacement per iteration is clamped to 30% of the polygon's equivalent radius to prevent topology-breaking jumps. Tobler (2004) reviewed the broader cartogram tradition and placed the diffusion method alongside Dorling's (1996) circular cartograms as the two most computationally tractable approaches.

Mathematical Formulation

Target area proportional to variable value. $$ A_i^{\text{target}} = v_i \cdot \frac{\sum_j A_j}{\sum_j v_j} $$ where $A_i$ is the current polygon area and $v_i$ is the area-representation field value. \tag{1}
Mass as equivalent-radius discrepancy. $$ r_i = \sqrt{\frac{A_i}{\pi}}, \qquad r_i^{\text{target}} = \sqrt{\frac{A_i^{\text{target}}}{\pi}}, \qquad m_i = r_i^{\text{target}} - r_i $$ Positive mass indicates the polygon needs to expand; negative mass indicates contraction. \tag{2}
Size error metric. $$ e_i = \frac{\max(A_i,\; A_i^{\text{target}})}{\min(A_i,\; A_i^{\text{target}})} $$ A value of 1.0 indicates exact area match. The algorithm halts when $\text{mean}(e_i) \leq \text{threshold}$. \tag{3}
Gaussian force kernel. $$ \mathbf{F}_{i \to j} = m_j \cdot r_j \cdot \gamma \cdot \exp\!\left(-\frac{d_{ij}^2}{2\sigma^2}\right) \cdot \hat{\mathbf{u}}_{ij} $$ where $d_{ij}$ is the distance between polygon centroids, $\sigma$ is the kernel bandwidth, $\gamma$ is a damping factor, and $\hat{\mathbf{u}}_{ij}$ is the unit vector from centroid $i$ to centroid $j$. \tag{4}
Displacement clamping. $$ \Delta_{\max} = r_i \cdot 0.3 $$ Per-iteration vertex displacement is capped at 30% of the polygon's equivalent radius to prevent topology tearing. \tag{5}

Parameters

ParameterTypeDefaultDescription
Input polygon layerPolygon vector—Enumeration areas to distort
Area-representation fieldNumeric field—The variable area should represent
Max iterationsInteger100Range 1–200; higher = more precise
Max avg error %Float5.0Range 0.1–100; halt criterion

Output Description

The output polygon layer contains the distorted geometries. Two metadata fields are appended: cartogram_error_pct (the final mean size error, expressed as a percentage of the target) and the algorithm log reports the number of iterations actually run. All original attributes are preserved.

Interpretation Guide

Interpretation A continuous-area cartogram replaces the "land area" mental model with a "variable weight" model: expanded regions indicate over-representation of the mapped variable relative to geographic area; shrunken regions indicate under- representation. The map remains topologically recognisable if the variable is at least loosely correlated with geographic area. For highly skewed variables (e.g., GDP where one city dominates), expect extreme distortion of the dominant unit. Pair the cartogram with a standard choropleth to help the reader calibrate their spatial expectations.

References

  1. Gastner, M. T. & Newman, M. E. J. Diffusion-Based Method for Producing Density-Equalizing Maps. Proceedings of the National Academy of Sciences, 101(20): 7499–7504, 2004. doi:10.1073/pnas.0400280101
  2. Tobler, W. R. Thirty Five Years of Computer Cartograms. Annals of the Association of American Geographers, 94(1): 58–73, 2004. doi:10.1111/j.1467-8306.2004.09401004.x
  3. Dorling, D. Area Cartograms: Their Use and Creation. Concepts and Techniques in Modern Geography, 59. University of East Anglia, 1996.
  4. Gusein-Zade, S. M. & Tikunov, V. S. A New Technique for Constructing Continuous Cartograms. Cartography and Geographic Information Systems, 20(3): 167–173, 1993.
  5. Dougenik, J. A., Chrisman, N. R., & Niemeyer, D. R. An Algorithm to Construct Continuous Area Cartograms. The Professional Geographer, 37(1): 75–81, 1985.
  6. Henriques, R., Bação, F., & Lobo, V. Carto-SOM: Cartogram Creation Using Self-Organizing Maps. International Journal of Geographical Information Science, 23(4): 483–511, 2009.
  7. Nusrat, S. & Kobourov, S. The State of the Art in Cartograms. Computer Graphics Forum, 35(3): 619–642, 2016.
graticule_grid Map Reference

6. Graticule / Reference Grid

Theoretical Background

A graticule is the network of meridians and parallels that provides the coordinate reference frame for a map. Its correct construction is foundational: the spacing, labelling, and visual hierarchy of grid lines directly affect the map reader's ability to locate features and estimate distances. Robinson et al. (1995) devoted a full chapter to the graticule as a map element, emphasising that the interval should be chosen so that the grid is visible but not dominant, typically producing 6–10 lines across the map extent.

The CartoLab implementation generates a vector line layer of meridians and parallels on "nice" intervals—round numbers drawn from the set $\{1, 2, 5\} \times 10^k$ that divide the coordinate space into visually manageable segments. When the user leaves an interval at 0, the algorithm computes the span of the extent in decimal degrees, divides by a target of approximately 8 lines, and normalises the result to the nearest nice number. Each line is tagged with its orientation (meridian/parallel), its constant coordinate value, and a formatted human-readable label. Lines are aligned to the global coordinate origin (multiples of the interval starting from $\lfloor x_{\min} / \Delta x \rfloor \cdot \Delta x$), ensuring that adjacent map sheets will share coincident grid lines.

Mathematical Formulation

Nice interval algorithm. $$ \text{nice}(s) = \min_{k \in \mathbb{Z}}\; \min_{b \in \{1,2,5\}}\; \left| b \times 10^{\lfloor \log_{10} s \rfloor + k} - s \right| $$ where $s = \text{span} / \text{target\_lines}$. The algorithm searches for the nearest round number in the geometric sequence $\{1, 2, 5\} \times 10^n$. \tag{1}
Aligned meridian positions. $$ x_k = \left\lceil \frac{x_{\min}}{\Delta x} \right\rceil \cdot \Delta x $$ for $k = 0, 1, \ldots$ while $x_k \leq x_{\max}$. \tag{2}
Parallel positions. $$ y_k = \left\lceil \frac{y_{\min}}{\Delta y} \right\rceil \cdot \Delta y $$ for $k = 0, 1, \ldots$ while $y_k \leq y_{\max}$. \tag{3}
Auto-interval computation. $$ \Delta x_{\text{auto}} = \operatorname{nice}\!\left( \frac{x_{\max} - x_{\min}}{8} \right), \qquad \Delta y_{\text{auto}} = \operatorname{nice}\!\left( \frac{y_{\max} - y_{\min}}{8} \right) $$ \tag{4}

Parameters

ParameterTypeDefaultDescription
Grid extentExtent (xmin, xmax, ymin, ymax)Canvas extentCRS-dependent; in degrees for geographic CRS
X intervalFloat0 (auto)Meridian spacing in CRS units; 0 = auto
Y intervalFloat0 (auto)Parallel spacing in CRS units; 0 = auto

Output Description

A line layer where each feature is a single meridian or parallel spanning the full extent. Attributes include orientation (string: “meridian” or “parallel”), coord (float, the constant coordinate value), and label (string, formatted with degree symbols and hemisphere suffix where appropriate). Lines are styled with a thin, low-opacity stroke suitable for map overlay.

Interpretation Guide

Interpretation The graticule is a reference element, not a data layer. Keep the stroke thin (0.15–0.3 pt) and colour neutral (grey or muted tone from the map palette). The auto-interval produces grids with roughly 8–10 lines across the map, which strikes a good balance for A4/letter layouts. For inset maps or overview maps, manually set larger intervals to avoid clutter. The grid should always be the bottom-most layer in the map stack so that data layers draw on top.

References

  1. Robinson, A. H., Morrison, J. L., Muehrcke, P. C., Kimerling, A. J., & Guptill, S. C. Elements of Cartography (6th ed.). Wiley, 1995. ISBN 978-0471555797.
  2. Snyder, J. P. Map Projections—A Working Manual. U.S. Geological Survey Professional Paper 1395, 1987.
  3. Canters, F. & Decleir, H. The World in Perspective: A Directory of World Map Projections. Wiley, 1989. ISBN 978-0471921479.
  4. Iliffe, J. & Lott, R. Datums and Map Projections for Remote Sensing, GIS, and Surveying (2nd ed.). Whittles Publishing, 2008. ISBN 978-1904445470.
  5. Krygier, J. & Wood, D. Making Maps: A Visual Guide to Map Design for GIS (3rd ed.). Guilford Press, 2016. ISBN 978-1462509980.
  6. Slocum, T. A. et al. Thematic Cartography and Geovisualization (3rd ed.). Pearson, 2009. ISBN 978-0132298346.
value_by_alpha Thematic Mapping

7. Value-by-Alpha (VbA) Map

Theoretical Background

The Value-by-Alpha (VbA) map, introduced by Roth, Woodruff, and Johnson (2010), addresses a persistent cartographic tension: how to encode two variables on a choropleth when one variable represents the magnitude of a phenomenon and the other represents the reliability or confidence in that measurement. Traditional bivariate maps treat both variables symmetrically, but reliability is not a symmetric partner to magnitude—it is a meta-property that should qualify, not compete with, the primary signal.

The VbA solution encodes the primary variable through colour (hue and saturation) and the reliability variable through opacity (alpha channel). High-reliability features are rendered fully opaque, drawing the reader's attention with saturated colour; low-reliability features are partially transparent, visually receding into the background. This perceptual hierarchy mirrors the information hierarchy: the reader should first notice the pattern of the primary variable and then adjust their confidence based on the visibility of each unit.

CartoLab implements VbA by classifying the reliability variable into $n$ bins using Fisher-Jenks natural breaks (to respect the distributional structure of reliability), then mapping class bins linearly to an alpha range specified by the user. The primary variable is styled independently (graduated colour); the VbA alpha field is added as a data-defined opacity override in the renderer. MacEachren (1995) and Bertin (1967/1983) provided the theoretical foundation for visual variables and their perceptual ordering, of which the CartoLab implementation directly applies the hierarchy of colour over transparency.

Mathematical Formulation

Alpha mapping from reliability bins. $$ \alpha_i = \alpha_{\min} + b_i \cdot \frac{\alpha_{\max} - \alpha_{\min}}{n_{\text{bins}} - 1} $$ where $b_i \in \{0, 1, \ldots, n_{\text{bins}}-1\}$ is the reliability bin index (0 = least reliable / lowest alpha) and $\alpha \in [0, 255]$ is the 8-bit alpha channel. \tag{1}
Fisher-Jenks classification of reliability. $$ \{t_0, t_1, \ldots, t_K\} = \operatorname{argmin} \sum_{k=1}^{K} \sum_{x_i \in [t_{k-1}, t_k)} (x_i - \bar{x}_k)^2 $$ where $t_0 = \min(x)$, $t_K = \max(x)$, and $K = n_{\text{bins}}$ is the number of alpha bins. \tag{2}
Opacity percentage (human-readable). $$ \alpha_i^{\text{pct}} = \frac{\alpha_i}{255} \times 100 $$ Rounded to one decimal place for the vba_alpha_pct output field. \tag{3}
Visual variable hierarchy (Bertin, 1967). The VbA map exploits the ordered visual variables: colour saturation $\succ$ transparency $\succ$ size $\succ$ orientation. By assigning the primary variable to saturation/hue (high in the perceptual hierarchy) and reliability to transparency (lower), the map respects the intended information priority. \tag{4}

Parameters

ParameterTypeDefaultDescription
Input layerPolygon vector—Must contain both primary and reliability fields
Primary variable (colour)Numeric field—Encoded through hue/saturation
Reliability variable (opacity)Numeric field—Higher = more opaque (typically)
Alpha minInteger40Range 0–255; opacity for the least reliable bin
Alpha maxInteger255Range 0–255; opacity for the most reliable bin

Output Description

The output polygon layer retains all original fields and adds vba_alpha (integer, 0–255, the computed alpha channel value) and vba_alpha_pct (float, 0–100, the same value as a human-readable percentage). These fields can be used as data-defined overrides for the opacity property in any QGIS renderer.

Interpretation Guide

Interpretation VbA maps require an explicit legend that explains both channels, ideally a two-part legend: a graduated colour ramp for the primary variable and a graduated alpha strip for reliability. Alpha min should not fall below 40 (0–255 scale) or features become effectively invisible. For print, test the minimum alpha on the target paper stock and press, as ink gain can alter the perceived transparency. The technique is most effective when reliability varies meaningfully across the map; if all units have similar reliability, the map reduces to an ordinary choropleth.

References

  1. Roth, R. E., Woodruff, A. W., & Johnson, Z. F. Value-by-Alpha Maps: An Alternative Technique to the Cartogram. The Cartographic Journal, 47(2): 130–140, 2010. doi:10.1179/000870409X12488753453372
  2. MacEachren, A. M. How Maps Work: Representation, Visualization, and Design. Guilford Press, 1995 (revised 2004). ISBN 978-1572300408.
  3. Bertin, J. Semiology of Graphics: Diagrams, Networks, Maps. University of Wisconsin Press, 1983. (Original French edition 1967.) ISBN 978-0299090603.
  4. Sun, H. & Li, Z. Effectiveness of Cartogram for the Representation of Spatial Data. The Cartographic Journal, 47(1): 12–21, 2010.
  5. Garlandini, S. & Fabrikant, S. I. Evaluating the Effectiveness and Efficiency of Visual Variables for Geographic Information Visualization. In: Spatial Information Theory (COSIT), LNCS 5756: 195–211. Springer, 2009.
  6. Brewer, C. A. & Pickle, L. Evaluation of Methods for Classifying Epidemiological Data on Choropleth Maps in a Series. Annals of the Association of American Geographers, 92(4): 662–681, 2002.
  7. Robinson, A. C. Visual Highlighting Methods for Geovisualization. Cartography and Geographic Information Science, 38(4): 373–383, 2011.
quick_style Quick Style

8. Quick Style

Theoretical Background

Map styling is the bridge between data and visual communication, yet the default single-symbol renderer in GIS software forces the analyst to manually configure classification, colour ramps, and legend formatting before any meaningful pattern is visible. Quick Style automates this pipeline with cartographically informed defaults, applying Harrower and Brewer's (2003) ColorBrewer palettes for categorical data and the perceptually uniform Viridis family for sequential numeric data. The goal is to produce a publication-ready styled layer in a single operation.

The algorithm detects the field type (numeric or string) and auto-selects a graduated or categorised renderer accordingly, with a force-mode override for cases where coded integers should be treated as categories. For graduated renderers, the user selects from three break methods: quantile (equal count per class), equal interval (equal data range per class), or geometric interval (optimal for skewed distributions, as described in Algorithm 1). The quantile method is the default for choropleth mapping per Brewer's (1994) recommendation for ordinal color schemes, as it ensures each colour class contains approximately the same number of observations.

The output includes a thin white outline (0.15 mm) that visually separates adjacent polygons without drawing attention to boundaries—a technique borrowed from print cartography where a hairline gap between fill areas improves figure-ground discrimination. The style summary string provides a compact report of the applied configuration for reproducibility.

Mathematical Formulation

Quantile breaks. Given sorted values $x_{(1)} \leq x_{(2)} \leq \cdots \leq x_{(n)}$, class boundaries are $$ b_k = x_{(\lceil k \cdot n / K \rceil)}, \qquad k = 1, \ldots, K-1 $$ with $b_0 = x_{(1)}$ and $b_K = x_{(n)}$. Each class contains $\lfloor n/K \rfloor$ or $\lceil n/K \rceil$ observations. \tag{1}
Equal-interval breaks. $$ b_k = x_{\min} + k \cdot \frac{x_{\max} - x_{\min}}{K}, \qquad k = 0, \ldots, K $$ \tag{2}
Geometric interval (as in Algorithm 1). $$ C_k = a \cdot r^k, \qquad (r^*, a^*) = \operatorname{argmin}_{r,a} \sum_{j=1}^{K} \left( f_j - \frac{n}{K} \right)^2 $$ \tag{3}
ColorBrewer palette assignment. For a categorical field with $m$ unique values, the first $m$ colours from the selected ColorBrewer qualitative palette are mapped 1:1 to the sorted unique values. For sequential data, $K$ colours are linearly interpolated from the selected sequential or diverging scheme. \tag{4}

Parameters

ParameterTypeDefaultDescription
Vector layerLoaded vector layer—Must be loaded in the current QGIS project
FieldAny field—Numeric triggers graduated; text triggers categorised
ModeEnumAutoAuto / Graduated / Categorised
ClassesInteger5Range 2–12
Break methodEnumQuantileQuantile / Equal Interval / Geometric Interval
PaletteEnumViridisColorBrewer or Viridis family palette
ReverseBooleanFalseInvert the colour ramp direction
Thin white outlineBooleanTrueAdds a 0.15 mm white stroke to polygon fills

Output Description

The output is an in-place styled layer: no new layer is created; the existing layer's renderer is replaced with the computed graduated or categorised renderer. The algorithm returns a summary string (e.g., “Graduated: pop_density, 5 quantile classes, Viridis”). The legend is automatically populated with formatted class labels.

Interpretation Guide

Interpretation Quick Style is a productivity tool, not a classification tool. For exploratory analysis, accept the auto defaults; for publication, use Quick Style to establish a baseline and then refine the renderer in the layer styling panel. The thin white outline is especially valuable for dark fills and when the map will be reproduced in greyscale, where adjacent polygons of similar value can merge.

References

  1. Harrower, M. & Brewer, C. A. ColorBrewer.org: An Online Tool for Selecting Colour Schemes for Maps. The Cartographic Journal, 40(1): 27–37, 2003. doi:10.1179/000870403235002042
  2. Brewer, C. A. Color Use Guidelines for Mapping and Visualization. In: MacEachren, A. M. & Taylor, D. R. F. (eds.) Visualization in Modern Cartography. Pergamon, 1994: 123–147.
  3. Brewer, C. A. Designing Better Maps: A Guide for GIS Users (2nd ed.). ESRI Press, 2015. ISBN 978-1589484405.
  4. Smith, R. M. Comparing Traditional Methods for Selecting Class Intervals on Choropleth Maps. The Professional Geographer, 38(1): 62–67, 1986.
  5. Monmonier, M. How to Lie with Maps (3rd ed.). University of Chicago Press, 2018. ISBN 978-0226435923.
  6. Kimerling, A. J., Buckley, A. R., Muehrcke, P. C., & Muehrcke, J. O. Map Use: Reading, Analysis, Interpretation (8th ed.). ESRI Press, 2016. ISBN 978-1589484429.
proportional_symbols Thematic Mapping

9. Proportional Symbols (Flannery)

Theoretical Background

Proportional symbol maps encode a quantitative variable through the size of a point symbol placed at each feature's centroid. The technique is particularly suited to count or magnitude data where the absolute value, rather than the density, is the quantity of interest. The central perceptual challenge is that map readers systematically underestimate the area of larger circles: Flannery (1971) demonstrated experimentally that the perceived size $S_p$ follows Stevens' (1957) power law, $S_p \propto S^{k}$ with $k \approx 0.8747$, rather than the geometric $S_p \propto S^{1.0}$ or $S_p \propto S^{0.5}$ (radius scaling). The common Flannery compensation exponent of $k = 0.5716$ adjusts the area scaling so that the reader's perceived ratio between symbols approximates the true data ratio.

The CartoLab implementation scales symbol size (in map millimetres) from the user-specified minimum and maximum, applying either Flannery-compensated scaling ($k = 0.5716$) or true area scaling ($k = 0.5$, i.e., radius $\propto \sqrt{v}$). The algorithm also computes suggested legend values at geometrically spaced intervals rounded to nice numbers (1, 2, 5 $\times 10^m$), providing the map author with a ready-made set of reference symbols. Cleveland and McGill (1984) ranked position along a common scale as the most accurate visual variable, but among area-based encodings, circle size is the most accurate and the most commonly used.

Mathematical Formulation

Size mapping with Flannery compensation. $$ s_i = s_{\min} + (s_{\max} - s_{\min}) \cdot \left( \frac{v_i}{v_{\max}} \right)^k $$ where $s_i$ is the symbol size in millimetres, $k = 0.5716$ (Flannery) or $k = 0.5$ (true area scale). When $k = 0.5$, the symbol area is linearly proportional to the data value. \tag{1}
Nice legend values (geometric spacing). Legend values $\{v_0, v_1, \ldots, v_{L-1}\}$ are chosen such that $$ v_j = \operatorname{nice}\!\left( v_{\max} \cdot g^{\,j - L + 1} \right), \qquad g = \left( \frac{v_{\max}}{v_{\min}} \right)^{1/(L-1)} $$ where $\operatorname{nice}(x)$ rounds $x$ to the nearest number of the form $d \times 10^m$ with $d \in \{1, 2, 5\}$. \tag{2}
Stevens' power law for perceived area. $$ \psi(S) = k \cdot S^{\,n}, \qquad n \approx 0.7 \text{ (for area perception)} $$ where $\psi$ is the perceived magnitude and $S$ is the physical stimulus magnitude. The Flannery exponent combines the area-perception exponent with the radius-to-area relationship. \tag{3}
The perceived ratio under Flannery compensation. $$ \frac{\psi(s_i)}{\psi(s_j)} \approx \frac{v_i}{v_j} $$ for any pair of symbols $(i, j)$, correcting the systematic underestimation of large circle areas relative to small ones. \tag{4}

Parameters

ParameterTypeDefaultDescription
Input layerAny vector—A point symbol is placed at each feature's centroid
Magnitude fieldNumeric field—Drives the symbol size; must be positive
Max symbol size (mm)Float10.0Symbol size for the maximum value
Min symbol size (mm)Float2.0Symbol size for the minimum value
Flannery compensationBooleanTrueApply $k = 0.5716$; false uses $k = 0.5$

Output Description

The output is a point layer with one point per input feature (at its centroid). New fields: psym_value (the original magnitude value, echoed for reference), psym_size (float, symbol diameter in millimetres). The layer is auto-styled with a data-defined marker size using psym_size. The algorithm log prints suggested legend values for manual legend construction.

Interpretation Guide

Interpretation Enable Flannery compensation for maps intended for a general audience, as it compensates for the well-documented tendency to underestimate large-circle areas. For scientific and professional audiences who may mentally calibrate against the legend, the true-area ($k = 0.5$) option is defensible. Always include at least three reference circles in the legend (minimum, median, maximum); the suggested legend values printed in the log provide geometrically spaced intermediate references. Overlapping symbols are a known limitation: for densely packed features, consider using a smaller maximum size or a cartogram approach instead.

References

  1. Flannery, J. J. The Relative Effectiveness of Some Common Graduated Point Symbols in the Presentation of Quantitative Data. The Canadian Cartographer, 8(2): 96–109, 1971.
  2. Stevens, S. S. On the Psychophysical Law. Psychological Review, 64(3): 153–181, 1957. doi:10.1037/h0046162
  3. Cleveland, W. S. & McGill, R. Graphical Perception: Theory, Experimentation, and Application to the Development of Graphical Methods. Journal of the American Statistical Association, 79(387): 531–554, 1984. doi:10.1080/01621459.1984.10478080
  4. Chang, K. Visualizing Geometric Uncertainty of Proportional Symbol Maps. Cartography and Geographic Information Science, 36(4): 349–365, 2009.
  5. Groop, R. E. & Cole, D. Overlapping Graduated Circles: Magnitude Estimation and Method of Portrayal. The Canadian Cartographer, 15(2): 114–122, 1978.
  6. Slocum, T. A. et al. Thematic Cartography and Geovisualization (3rd ed.). Pearson, 2009. ISBN 978-0132298346.
  7. Tanimura, S., Kuroiwa, C., & Mizota, T. Proportional Symbol Mapping in R. Journal of Statistical Software, 15(5): 1–7, 2006.
ridge_map Thematic Mapping

10. Ridge Map (Joy Division Style)

Theoretical Background

The ridge map, sometimes called a Joy Division map after Peter Saville's iconic album cover design for Unknown Pleasures (which itself adapted a plot of radio pulsar CP 1919), transforms a raster surface into a set of vertically deformed scanlines. Each scanline is a horizontal transect across the raster; pixel values displace the line vertically, creating an overlapping sequence of wave-like profiles. The resulting graphic encodes elevation, density, or any continuous raster variable as a rhythmic spatial texture that is immediately legible as topography.

The technique owes its intellectual lineage to Tufte's (1983) principle of data-ink maximisation and Imhof's (1982) cartographic relief presentation, both of which argue that surface representation should reveal the underlying data structure with minimal decorative intervention. The ridge map's aesthetic of overlapping transparent lines produces a figure-ground effect where crests (high pixel values) visually advance and troughs (low values) recede, mimicking shaded relief despite using only vector linework.

CartoLab implements ridge mapping by sampling the input raster along equally spaced horizontal scanlines, deforming each $y$-coordinate by the pixel value multiplied by a user-specified vertical exaggeration. An optional Laplacian smoothing pass (averaging each vertex with its two neighbours) reduces high-frequency noise from the raster. The scanlines are clipped to an optional user-supplied polygon extent. The output is a vector line layer suitable for rendering with dark-background and transparency for the overlapping ridge effect.

Mathematical Formulation

Scanline deformation. $$ y'(x) = y_{\text{baseline}} + z(x, y_{\text{baseline}}) \cdot \lambda $$ where $z(x, y)$ is the raster pixel value at column $x$, row $y$, and $\lambda$ is the vertical exaggeration factor. $y_{\text{baseline}}$ is the un-deformed scanline $y$-coordinate. \tag{1}
Laplacian smoothing (one pass). $$ p_i' = \frac{p_{i-1} + p_i + p_{i+1}}{3} $$ Applied along the scanline vertex sequence. Repeated $m$ times for $m$ smoothing passes (0–20). \tag{2}
Scanline spacing. $$ \Delta y_{\text{line}} = \frac{y_{\max} - y_{\min}}{n_{\text{lines}} - 1} $$ where $n_{\text{lines}}$ is the user-specified number of scanlines spanning the raster's vertical extent. \tag{3}
Raster value sampling (bilinear interpolation). $$ z(x, y) = (1 - u)(1 - w)\,z_{00} + u(1 - w)\,z_{10} + (1 - u)w\,z_{01} + u w\,z_{11} $$ where $u, w \in [0, 1]$ are the fractional column and row offsets within the enclosing four-pixel neighbourhood. \tag{4}

Parameters

ParameterTypeDefaultDescription
Input raster (single-band)Raster layer—Single-band; multi-band uses the first band
Number of scanlinesInteger100Range 5–500; more = denser ridge pattern
Vertical exaggerationFloat1.0Multiplier on pixel values in map units
Line spacingFloat00 = auto (equal spacing); manual override in map units
Smoothing passesInteger1Range 0–20; 0 = raw raster values
Clip extentPolygon (optional)—Restricts output to a specific region

Output Description

A LineString vector layer with line_id (integer, scanline index from top to bottom). Each line contains the deformed vertex sequence. The recommended styling uses a dark background (#132025 or similar), a light stroke (white or pale colour), low opacity (10–40%), and no fill. The line_id field enables per-line styling, e.g., data-defined colour ramps.

Interpretation Guide

Interpretation Ridge maps work best with continuous raster surfaces: DEMs, density rasters, and kernel-density outputs. The visual effect relies on overlapping transparent lines: on a dark background, light lines accumulate brightness in valleys (where many lines converge at similar $y'$) and remain dim on ridges (where lines diverge). Start with 100 scanlines and 10% stroke opacity; increase scanlines for smoother surfaces and decrease for a sparser, more "drawn" aesthetic. Smoothing above 5 passes begins to suppress real topographic features; use sparingly.

References

  1. Tufte, E. R. The Visual Display of Quantitative Information. Graphics Press, 1983 (2nd ed. 2001). ISBN 978-1930824133.
  2. Imhof, E. Cartographic Relief Presentation. De Gruyter, 1982 (English ed. 2007, edited by H. J. Steward). ISBN 978-1589484580.
  3. Wood, D. The Power of Maps. Guilford Press, 1992. ISBN 978-0898624922.
  4. Lobeck, A. K. Block Diagrams and Other Graphic Methods Used in Geology and Geography. Emerson-Trussell, 1924 (2nd ed. 1958).
  5. Tanaka, K. The Relief Contour Method of Representing Topography on Maps. Geographical Review, 40(3): 444–456, 1950.
  6. Kennelly, P. J. & Kimerling, A. J. Desktop Hachure Maps from Digital Elevation Models. Cartographic Perspectives, 37: 78–81, 2000.
  7. Brewer, C. A. Designing Better Maps: A Guide for GIS Users (2nd ed.). ESRI Press, 2015. ISBN 978-1589484405.
building_25d_style 2.5D Styling

11. Building 2.5D Style

Theoretical Background

The 2.5D representation bridges the gap between two-dimensional planimetric maps and full three-dimensional scenes. By extruding polygon footprints vertically and rendering them with an isometric or pseudo-3D projection, the technique conveys building height, massing, and urban form while retaining the precision and measurability of a 2D map. QGIS' native 2.5D renderer, introduced in version 3.x, applies an isometric affine transformation to extruded geometry and renders walls and roofs with user-configurable materials and lighting.

The CartoLab implementation offers two modes: the native QGIS 2.5D renderer (configuring angle, height scale, shadow, wall colour, and roof colour via a live preset system) and the per-floor colour-band mode. The per-floor mode divides the building height into floor-height increments and applies a repeating colour palette to each floor band, producing a legend-friendly single-symbol renderer where each floor is a separate symbol layer. This approach, inspired by Jenny et al.'s (2010) design principles for Swiss-style rock drawing, uses systematic colour coding to reveal floor patterns across an urban landscape.

The visual presets (Warm Civic, Cool Slate, Limestone & Teal, Night Copper, Urban Core Gold, Suburban Pastoral, Cyberpunk Night, Nordic Minimalist, Emerald Metropolis, and Terracotta Mediterranean) apply pre-tuned combinations of roof colour, wall colour, shadow strength, and material style. The height can be supplied in metres (direct field value) or in floor count (field value $\times$ floor height in metres). An optional height clamp caps extreme outliers, and stepped extrusion toggles between smooth (continuous extrusion) and stepped (floor-level increments) wall profiles.

Mathematical Formulation

Height computation. $$ h_i = \begin{cases} v_i, & \text{metres mode} \\[4pt] v_i \cdot h_{\text{floor}}, & \text{floors mode} \end{cases} $$ where $v_i$ is the height field value for feature $i$ and $h_{\text{floor}}$ is the floor height in metres. \tag{1}
2.5D isometric projection. $$ \mathbf{T} = \begin{bmatrix} \cos\theta & -\cos\theta & 0 \\ -\sin\theta \cdot \sin\phi & -\sin\theta \cdot \sin\phi & \cos\phi \end{bmatrix} $$ where $\theta$ is the rotation angle (user-specified, typically 45° or 30°), and $\phi$ is the pitch (typically 45° for isometric). \tag{2}
Floor band height assignment. $$ n_i^{\text{floors}} = \left\lceil \frac{h_i}{h_{\text{floor}}} \right\rceil $$ Each floor $k \in \{1, \ldots, n_i^{\text{floors}}\}$ receives colour $\text{palette}[k \bmod |\text{palette}|]$ from repeating palettes (Civic Spectrum, Planning Bands, Soft Atlas, Turbo Height, or Viridis Height). \tag{3}
Height clamp. $$ h_i' = \min(h_i, h_{\max}) $$ Applied after height mode computation; features exceeding $h_{\max}$ are capped. \tag{4}
Shadow offset (pseudo-lighting). A parallel shadow polygon is displaced by $\Delta x = s \cdot \cos\theta_s$, $\Delta y = s \cdot \sin\theta_s$, where $s$ is the shadow distance (typically proportional to building height) and $\theta_s$ is the light azimuth. \tag{5}

Parameters

ParameterTypeDefaultDescription
Polygon layerLoaded polygon vector—Building footprint layer
Height fieldNumeric field—Height in metres or floor count
Height modeEnumMetresMetres (direct value) or Floors (count * floor height)
Floor heightFloat3.0Metres per floor; used in floors mode and floor-band renderer
Renderer styleEnum2.5DNative 2.5D renderer or Per-floor colour bands
Floor paletteEnumCivic SpectrumCivic Spectrum, Planning Bands, Soft Atlas, Turbo Height, Viridis Height
Max floorsInteger30Maximum number of symbol layers generated
Visual presetEnumWarm CivicWarm Civic, Cool Slate, Limestone & Teal, Night Copper, Urban Gold, Suburban Pastoral, Cyberpunk, Nordic, Emerald Metropolis, Terracotta Mediterranean
AngleFloat110.0Viewing azimuth in degrees
Height scaleFloat1.0Vertical exaggeration multiplier
Max height clampFloat0 (off)Cap extreme heights; 0 = no clamp
SteppedBooleanFalseStepped extrusion at floor intervals
ShadowBooleanTrueCast shadow from a 315° light source
Wall shadingBooleanTrueDarken walls facing away from the light

Output Description

The output is an in-place styled layer: the existing layer's renderer is replaced. A summary string describes the applied configuration (e.g., “2.5D: Warm Civic preset, 110° angle, height_metres, shadow on”). For the floor-band mode, a separate symbol layer is created for each floor level up to Max Floors, each with a data-defined height filter.

Interpretation Guide

Interpretation The 2.5D view is a qualitative representation; do not use it for precise height measurement. The floor-band mode excels at visualising building-height policy zones (e.g., where 2–4 floor buildings abut 10+ floor towers) and works well in print layouts. Choose a preset that contrasts with the basemap: Emerald Metropolis for environmental and eco-district master plans, Terracotta Mediterranean for historic old town quarters, and Cool Slate for general architectural master plans.

References

  1. QGIS Development Team. QGIS 2.5D Renderer Documentation. QGIS User Manual, 2024. qgis.org/docs
  2. Jenny, B., Gilgen, J., Geisthövel, R., Marston, B. E., & Hurni, L. Design Principles for Swiss-Style Rock Drawing. The Cartographic Journal, 47(4): 323–332, 2010.
  3. Kraak, M. J. & Ormeling, F. Cartography: Visualization of Geospatial Data (3rd ed.). Guilford Press, 2011. ISBN 978-1609181949.
  4. Haeberling, C., Bär, H., & Hurni, L. Proposed Cartographic Design Principles for 3D Maps. Cartography and Geographic Information Science, 35(4): 309–324, 2008.
  5. Döllner, J. & Buchholz, H. Continuous Level-of-Detail Modeling of Buildings in 3D City Models. Proceedings of ACM GIS, 2005.
  6. Jobst, M. & Döllner, J. Better Perception of 3D-Spatial Relations by View-Dependent Mixed-Mode Rendering of 2D and 3D Maps. International Journal of Geographical Information Science, 22(10): 1085–1103, 2008.
  7. Kwan, M. P. & Lee, J. Geovisualization of Human Activity Patterns Using 3D GIS: A Time-Geographic Approach. In: Goodchild, M. F. & Janelle, D. G. (eds.) Spatially Integrated Social Science. Oxford University Press, 2004.
label_points Labeling

12. Visual-Center Label Points

Theoretical Background

The visual centre of a polygon—the point from which a label appears most naturally to belong to the enclosing area—is not the centroid. Centroids can fall outside their polygon for concave shapes (e.g., C-shaped administrative units, crescent-shaped lakes, doughnut polygons with interior holes). The cartographic requirement is for a point that is guaranteed to lie inside the polygon and is as far as possible from any edge, maximising the available space for a label.

The polylabel algorithm, developed by Mapbox (2016), solves this problem efficiently using a quad-tree search with a priority queue. The algorithm maintains a max-heap ordered by each cell's potential: the distance from the cell centre to the nearest polygon edge plus the cell's half-diagonal (the maximum possible improvement from subdividing that cell). At each step, the cell with the highest potential is subdivided into four children, and the candidate point is updated whenever a cell centre has a larger signed distance than the current best. The search terminates when the best cell's potential falls within the specified precision of the best known distance.

The signed distance uses ray-casting to determine whether a point is inside the polygon; points outside receive a negative distance. This approach generalises to polygons with holes, multi-part geometries, and complex concavities. Garcia-Castellanos and Lombardo (2007) provided the theoretical analogue for the "pole of inaccessibility" on Earth's land surface; the polylabel algorithm adapts the quad-tree search for the vector case.

Mathematical Formulation

Signed distance to polygon boundary. $$ d(p) = \begin{cases} +\min_{e \in E} \operatorname{dist}(p, e), & \text{if } p \text{ is inside the polygon} \\[4pt] -\min_{e \in E} \operatorname{dist}(p, e), & \text{if } p \text{ is outside} \end{cases} $$ where $E$ is the set of boundary edges and inside/outside is determined by ray-casting (even-odd rule). \tag{1}
Point-to-segment distance squared. For edge $e$ from $a$ to $b$, the squared distance to point $p$ is: $$ d^2(p, e) = \begin{cases} |p - a|^2, & \text{if } (p-a)\cdot(b-a) \leq 0 \\[4pt] |p - b|^2, & \text{if } (p-b)\cdot(a-b) \leq 0 \\[4pt] \dfrac{|(p-a) \times (b-a)|^2}{|b-a|^2}, & \text{otherwise} \end{cases} $$ where $\times$ denotes the 2D cross product. \tag{2}
Cell potential (quad-tree priority). $$ \operatorname{potential}(C) = d(\text{centre}_C) + h_C \cdot \sqrt{2} $$ where $h_C$ is the half-side length of cell $C$. This is the maximum possible signed distance achievable by any point within the cell. \tag{3}
Quad-tree termination condition. The search halts when $$ \max_{C \in Q} \operatorname{potential}(C) - d(p^*) \leq \varepsilon $$ where $p^*$ is the current best point, $Q$ is the priority queue, and $\varepsilon$ is the precision parameter. \tag{4}
Auto-precision scaling. $$ \varepsilon = \frac{\max(w, h)}{100} $$ where $w$ and $h$ are the width and height of the polygon's bounding box. \tag{5}
Principal Inertia Orientation Angle ($\text{lbl\_angle}$). $$ M_{xx} = \frac{1}{m}\sum_{i=1}^m (x_i - \bar{x})^2, \quad M_{yy} = \frac{1}{m}\sum_{i=1}^m (y_i - \bar{y})^2, \quad M_{xy} = \frac{1}{m}\sum_{i=1}^m (x_i - \bar{x})(y_i - \bar{y}) $$ $$ \theta_{\text{orient}} = \frac{1}{2}\operatorname{atan2}\bigl(2 M_{xy}, \, M_{xx} - M_{yy}\bigr) \times \frac{180^\circ}{\pi} \in [-90^\circ, +90^\circ] $$ Determines the dominant longitudinal elongation axis of parcels, river corridors, and elongated polygons for aligned text placement. \tag{6}

Parameters

ParameterTypeDefaultDescription
Input polygon layerPolygon vector—Polygons for which to compute label points
PrecisionFloat0 (auto)Search precision in map units; 0 auto-scales from bbox/100

Output Description

A point layer with one point per input polygon. Each point carries:

These fields enable data-defined label sizing and automatic shape-aligned rotation.

Interpretation Guide

Interpretation The label point's lbl_dist is the key quality metric: a value less than approximately 2 mm in map units at the target scale suggests that the polygon is too narrow for a readable label and should be labelled with a leader line or left unlabelled. The lbl_angle field provides the exact angle to align map labels with elongated parcels, building wings, or river segments.

References

  1. Mapbox. Polylabel: A Fast Algorithm for Finding the Pole of Inaccessibility of a Polygon. GitHub/mapbox/polylabel, 2016.
  2. Garcia-Castellanos, D. & Lombardo, U. Poles of Inaccessibility: A Calculation Algorithm for the Remotest Places on Earth. Scottish Geographical Journal, 123(3): 227–233, 2007. doi:10.1080/14702540801897856
  3. de Berg, M., Cheong, O., van Kreveld, M., & Overmars, M. Computational Geometry: Algorithms and Applications (3rd ed.). Springer, 2008. ISBN 978-3540779735.
  4. Edmondson, S., Christensen, J., Marks, J., & Shieber, S. A General Cartographic Labelling Algorithm. Cartographica, 33(4): 13–23, 1996.
normalize_field Data Preparation

13. Choropleth Normalization & Rates

Theoretical Background

Raw counts are nearly always the wrong variable for a choropleth map. Plotting total population, total crime incidents, or total disease cases on a map of unequal-area enumeration units produces a map of the size of the units, not the intensity of the phenomenon. Normalisation transforms a raw magnitude into a rate, density, index, or relative score that is comparable across spatial units. Monmonier (1991) devoted an entire chapter to the dangers of mapping un-normalised counts, and the American Cartographic Association's (Dent et al., 2009) guidelines consider normalisation a prerequisite for any choropleth intended for comparison.

CartoLab implements six normalisation methods spanning the principal use cases: (1) Rate divides a numerator by a denominator and scales the result, producing a density or proportion; (2) Z-score standardises to zero mean and unit variance, suitable for statistical comparison; (3) Robust z uses median and MAD (median absolute deviation) instead of mean and standard deviation, resisting the influence of outliers; (4) Min-max scales to the range $[0, 1]$, useful for input to multi-criteria analysis; (5) Percentile rank maps each value to its empirical cumulative probability $\times 100$, providing a purely ordinal scale; (6) $\log_{10}$ compresses the dynamic range of power-law or log-normal data, with an automatic shift to handle zero and negative values.

The choice of normalisation method is as consequential as the choice of classification. Krygier and Wood (2016) advise matching the normalisation to the analytical question: rates for comparisons of intensity, z-scores for identifying outliers, percentile ranks for ordinal rankings, and log transforms for visualising data spanning multiple orders of magnitude. The CartoLab implementation auto-graduates the output with the Viridis colour ramp and records the method in the output field for transparency.

Mathematical Formulation

1. Rate normalisation. $$ n_i = \frac{x_i}{d_i} \times S $$ where $d_i$ is the denominator (e.g., area, population) and $S$ is the rate scale (e.g., 100 for per cent, 1000 for per mille, 100000 for per 100k). \tag{1}
2. Z-score standardisation. $$ z_i = \frac{x_i - \mu}{\sigma}, \qquad \mu = \frac{1}{n}\sum_{i} x_i, \qquad \sigma = \sqrt{\frac{1}{n}\sum_{i} (x_i - \mu)^2} $$ \tag{2}
3. Robust z-score (median/MAD). $$ z_i^{\text{rob}} = \frac{x_i - \operatorname{median}(x)}{1.4826 \cdot \operatorname{MAD}}, \qquad \operatorname{MAD} = \operatorname{median}(|x_i - \operatorname{median}(x)|) $$ The constant $1.4826$ makes MAD a consistent estimator of $\sigma$ for normally distributed data. \tag{3}
4. Min-max normalisation. $$ n_i = \frac{x_i - x_{\min}}{x_{\max} - x_{\min}} \in [0, 1] $$ \tag{4}
5. Percentile rank. $$ p_i = 100 \times \frac{\operatorname{count}(x_j < x_i) + 0.5 \cdot \operatorname{count}(x_j = x_i)}{n} $$ Ties are handled by mid-rank assignment, producing values in $[0, 100]$. \tag{5}
6. Log-10 transform with shift. $$ n_i = \log_{10}(x_i + \delta), \qquad \delta = \begin{cases} 0, & \text{if } x_{\min} > 0 \\[4pt] |x_{\min}| + 1, & \text{if } x_{\min} \leq 0 \end{cases} $$ \tag{6}
7. Location Quotient (LQ / Specialisation Index). $$ LQ_i = \frac{e_i / E_i}{\sum_{k=1}^n e_k / \sum_{k=1}^n E_k} $$ where $e_i$ is local subsector activity and $E_i$ is total local employment or population. Values $> 1.0$ indicate regional concentration / competitive specialization. \tag{7}
8. Winsorized Min-Max Normalisation (Outlier-Clamped). $$ x_i^{\text{win}} = \max\bigl(P_5, \min(x_i, P_{95})\bigr), \qquad n_i = \frac{x_i^{\text{win}} - P_5}{P_{95} - P_5} \in [0, 1] $$ Clamps extreme $5^{\text{th}}$ and $95^{\text{th}}$ percentiles to prevent statistical blowout of map color gradients. \tag{8}
9. Decile Rank Normalisation ($1-10$). $$ D_i = \min\left(10, \left\lfloor 1 + \frac{10 \times \operatorname{rank}(x_i)}{n} \right\rfloor\right) \in \{1, 2, \dots, 10\} $$ Divides observations into 10 equal-frequency classes for socioeconomic ranking and indexing. \tag{9}

Parameters

ParameterTypeDefaultDescription
Input layerVector (any geometry)—Layer containing the field to normalise
Value fieldNumeric field—Raw count or magnitude to transform
MethodEnumRate0 = Rate, 1 = Z-score, 2 = Robust z, 3 = Min-max, 4 = Percentile rank, 5 = Log10, 6 = Location Quotient, 7 = Winsorized Min-Max, 8 = Decile Rank
Denominator fieldNumeric field—Required for Rate and Location Quotient methods; ignored for others
Rate scaleFloat100Multiplier applied to the ratio; e.g., 100 = %

Output Description

The output layer is a copy of the input with two new fields: norm_value (float, the transformed value) and norm_method (string, the method name applied, e.g., “location_quotient”). The layer is auto-graduated with the Viridis colour ramp. For the Rate and Location Quotient methods, denominator metrics are recorded in the processing log for scientific provenance.

Interpretation Guide

Interpretation Never map a raw count without a normalisation step. Rate is the standard default for choropleth mapping. Use Location Quotient ($LQ > 1$) to identify industrial clusters, creative clusters, and regional economic specialisation. Use Winsorized Min-Max when severe outliers compress the visual range of the remaining 90% of features. Use Decile Rank ($1-10$) for public policy and socio-economic deprivation indices.

References

  1. Monmonier, M. How to Lie with Maps (3rd ed.). University of Chicago Press, 2018. ISBN 978-0226435923.
  2. Dent, B. D., Torguson, J. S., & Hodler, T. W. Cartography: Thematic Map Design (6th ed.). McGraw-Hill, 2009. ISBN 978-0072943825.
  3. Isard, W. Methods of Regional Analysis: An Introduction to Regional Science. MIT Press, 1960.
  4. Tukey, J. W. Exploratory Data Analysis. Addison-Wesley, 1977.
  5. Rousseeuw, P. J. & Croux, C. Alternatives to the Median Absolute Deviation. Journal of the American Statistical Association, 88(424): 1273–1283, 1993.
bivariate_studio Interactive Studio

14. Interactive Bivariate Studio & Color Matrix Widget

Theoretical Background

The Interactive Bivariate Studio provides real-time exploration of 2D spatial relationships directly inside QGIS. By combining two continuous variables into an $N \times N$ matrix ($3\times3$ or $4\times4$), it reveals spatial correlations, spatial mismatches, and trade-offs (such as Social Vulnerability vs. Flood Hazard, or Public Transit Accessibility vs. Land Value).

Color Matrix Interpolation

Bilinear RGB blending in unit square $(u, v) \in [0, 1]^2$. $$ C(u, v) = (1-u)(1-v)C_{LL} + u(1-v)C_{HL} + (1-u)vC_{LH} + uvC_{HH} $$ where $C_{LL}$ is low-low (neutral grey), $C_{HL}$ is high-variable-1, $C_{LH}$ is high-variable-2, and $C_{HH}$ is high-high (convergent chromatic accent). \tag{1}

Curated Palette Themes

Preset$C_{LL}$$C_{LH}$$C_{HL}$$C_{HH}$Cartographic Use Case
Classic Brown-Teal#e8e8e8#5ab4ac#d8b365#8c510aGeneral socio-economic & environmental mapping
Pink-Blue (Robin)#e8e8e8#64acbe#c85a5a#574249Urban vulnerability & demographic shifts
Purple-Green#e8e8e8#7b3294#008837#2c103bLand use & ecological suitability
Warm Orange-Cyan#e8e8e8#2b83ba#fdae61#d7191cHeat island vs. green canopy coverage
Sunset Purple-Gold#e8e8e8#7b2cbf#f4a261#e76f51Housing affordability vs. employment access
print_layout_suite Layout Studio

15. Print Layout Studio Suite

Theoretical Background

The Print Layout Studio automates high-end cartographic layout production, eliminating tedious manual drafting while enforcing strict ISO page geometry, Swiss typography hierarchies, optical margins, and harmonious spatial balance.

Key Capabilities

color_accessibility Accessibility

16. Color Accessibility & CVD Simulation Engine

Theoretical Background

Ensures maps comply with WCAG 2.1 AAA accessibility standards and remain legible for individuals with Color Vision Deficiencies (CVD / Daltonism: Deuteranopia, Protanopia, Tritanopia).

Relative Luminance (sRGB). $$ L = 0.2126 R' + 0.7152 G' + 0.0722 B' $$ $$ C' = \begin{cases} \frac{C}{12.92}, & C \leq 0.04045 \\[4pt] \left(\frac{C + 0.055}{1.055}\right)^{2.4}, & C > 0.04045 \end{cases} \quad \text{for } C \in \{R, G, B\} / 255 $$ \tag{1}
WCAG Contrast Ratio. $$ CR = \frac{L_1 + 0.05}{L_2 + 0.05} \ge 7.0 \quad (\text{WCAG AAA for normal text}) $$ \tag{2}
sun_lighting Solar Lighting

17. Solar Position & Sun Angle Calculator

Theoretical Background

Computes real astronomical solar altitude and azimuth based on geographic latitude ($\phi$), day of year ($n$), and solar hour ($H$) to generate accurate 2.5D building shadows and daylight analysis.

Solar Declination & Altitude Equations. $$ \delta = 23.45^\circ \sin\left(\frac{360^\circ}{365}(284 + n)\right) $$ $$ \sin \alpha = \sin \phi \sin \delta + \cos \phi \cos \delta \cos H $$ $$ \cos \gamma = \frac{\sin \alpha \sin \phi - \sin \delta}{\cos \alpha \cos \phi} $$ where $\alpha$ is the solar altitude angle above the horizon and $\gamma$ is the solar azimuth angle. \tag{1}
scientific_palettes Cartographic Palettes

18. Perceptual Scientific & Topographic Color Ramps Library

Curated Scientific Ramps

PaletteTypeColorblind SafePrimary Use Case
TurboSequential Multi-HuePartialHigh-contrast continuous density & urban heat islands
MakoSequential MonotonicYesWater depth, marine morphology & wetlands
RocketSequential HotYesFire risk, building density & traffic congestion
EarthTopographic DivergingYesElevation gradients, land-use zoning & soil moisture
BathymetryMarine SequentialYesCoastal bathymetry & flood inundation levels
IceFireCyclic / DivergingYesTemperature anomalies, thermal stress & solar radiation