|
| 1 | +#' anyplot.ai |
| 2 | +#' map-projections: World Map with Different Projections |
| 3 | +#' Library: ggplot2 3.5.1 | R 4.4.1 |
| 4 | +#' Quality: 85/100 | Created: 2026-05-23 |
| 5 | + |
| 6 | +library(ggplot2) |
| 7 | +library(maps) |
| 8 | +library(ragg) |
| 9 | + |
| 10 | +# Theme tokens |
| 11 | +THEME <- Sys.getenv("ANYPLOT_THEME", "light") |
| 12 | +PAGE_BG <- if (THEME == "light") "#FAF8F1" else "#1A1A17" |
| 13 | +INK <- if (THEME == "light") "#1A1A17" else "#F0EFE8" |
| 14 | +INK_SOFT <- if (THEME == "light") "#4A4A44" else "#B8B7B0" |
| 15 | +OCEAN_BG <- if (THEME == "light") "#C4DCF0" else "#152030" |
| 16 | +LAND_BG <- if (THEME == "light") "#D4C9A8" else "#2A3020" |
| 17 | + |
| 18 | +ANYPLOT_PALETTE <- c( |
| 19 | + "#009E73", "#9418DB", "#B71D27", "#16B8F3", |
| 20 | + "#99B314", "#D359A7", "#BA843E" |
| 21 | +) |
| 22 | + |
| 23 | +# --- Mollweide projection math ----------------------------------------------- |
| 24 | +# Solve 2*theta + sin(2*theta) = pi*sin(lat) via Newton-Raphson |
| 25 | +solve_theta <- function(phi_vec) { |
| 26 | + sapply(phi_vec, function(p) { |
| 27 | + if (is.na(p)) return(NA_real_) |
| 28 | + if (abs(p) >= pi / 2 - 1e-9) return(sign(p) * pi / 2) |
| 29 | + t <- p |
| 30 | + for (i in seq_len(60)) { |
| 31 | + delta <- (2 * t + sin(2 * t) - pi * sin(p)) / (2 + 2 * cos(2 * t)) |
| 32 | + t <- t - delta |
| 33 | + if (abs(delta) < 1e-11) break |
| 34 | + } |
| 35 | + t |
| 36 | + }) |
| 37 | +} |
| 38 | + |
| 39 | +mollweide <- function(lon_deg, lat_deg) { |
| 40 | + lon <- lon_deg * pi / 180 |
| 41 | + lat <- lat_deg * pi / 180 |
| 42 | + theta <- solve_theta(lat) |
| 43 | + data.frame( |
| 44 | + x = (2 * sqrt(2) / pi) * lon * cos(theta), |
| 45 | + y = sqrt(2) * sin(theta) |
| 46 | + ) |
| 47 | +} |
| 48 | + |
| 49 | +# --- World country boundaries (projected) ------------------------------------ |
| 50 | +world <- map_data("world") |
| 51 | +world_proj <- mollweide(world$long, world$lat) |
| 52 | +world$x <- world_proj$x |
| 53 | +world$y <- world_proj$y |
| 54 | + |
| 55 | +# --- Projection boundary (Mollweide ellipse: a=2√2, b=√2) ------------------- |
| 56 | +t_ell <- seq(0, 2 * pi, length.out = 721) |
| 57 | +boundary <- data.frame( |
| 58 | + x = 2 * sqrt(2) * cos(t_ell), |
| 59 | + y = sqrt(2) * sin(t_ell) |
| 60 | +) |
| 61 | + |
| 62 | +# --- Graticule (lat/lon grid at 30° intervals) -------------------------------- |
| 63 | +lat_dense <- seq(-90, 90, length.out = 541) |
| 64 | +lon_dense <- seq(-180, 180, length.out = 1081) |
| 65 | + |
| 66 | +meridians <- do.call(rbind, lapply(seq(-180, 180, by = 30), function(lon0) { |
| 67 | + pts <- mollweide(rep(lon0, length(lat_dense)), lat_dense) |
| 68 | + data.frame(x = pts$x, y = pts$y, group = paste0("mer_", lon0)) |
| 69 | +})) |
| 70 | + |
| 71 | +parallels <- do.call(rbind, lapply(seq(-90, 90, by = 30), function(lat0) { |
| 72 | + pts <- mollweide(lon_dense, rep(lat0, length(lon_dense))) |
| 73 | + data.frame(x = pts$x, y = pts$y, group = paste0("par_", lat0)) |
| 74 | +})) |
| 75 | + |
| 76 | +graticule <- rbind(meridians, parallels) |
| 77 | + |
| 78 | +# --- Tissot indicatrices (small geodesic circles on the sphere) --------------- |
| 79 | +# Radius in degrees; longitude offset corrected for latitude (cos projection) |
| 80 | +r_deg <- 5 |
| 81 | +t_circle <- seq(0, 2 * pi, length.out = 73) # 73 pts → closed polygon |
| 82 | + |
| 83 | +lat_centers <- seq(-60, 60, by = 30) # 5 latitude bands |
| 84 | +lon_centers <- seq(-150, 150, by = 60) # 6 longitude positions |
| 85 | + |
| 86 | +tissot <- do.call(rbind, lapply(lat_centers, function(lat0) { |
| 87 | + cos_lat <- max(cos(lat0 * pi / 180), 0.08) |
| 88 | + do.call(rbind, lapply(lon_centers, function(lon0) { |
| 89 | + lon_c <- lon0 + r_deg * cos(t_circle) / cos_lat |
| 90 | + lat_c <- pmax(pmin(lat0 + r_deg * sin(t_circle), 89.9), -89.9) |
| 91 | + pts <- mollweide(lon_c, lat_c) |
| 92 | + data.frame( |
| 93 | + x = pts$x, |
| 94 | + y = pts$y, |
| 95 | + group = paste0("t_", lon0, "_", lat0) |
| 96 | + ) |
| 97 | + })) |
| 98 | +})) |
| 99 | + |
| 100 | +# --- Equator and prime meridian highlights ----------------------------------- |
| 101 | +equator <- mollweide(lon_dense, rep(0, length(lon_dense))) |
| 102 | +equator$group <- "equator" |
| 103 | + |
| 104 | +prime_merid <- mollweide(rep(0, length(lat_dense)), lat_dense) |
| 105 | +prime_merid$group <- "prime" |
| 106 | + |
| 107 | +# --- Plot -------------------------------------------------------------------- |
| 108 | +p <- ggplot() + |
| 109 | + # Ocean fill |
| 110 | + geom_polygon( |
| 111 | + data = boundary, aes(x = x, y = y), |
| 112 | + fill = OCEAN_BG, color = NA |
| 113 | + ) + |
| 114 | + # Land masses (country polygons) |
| 115 | + geom_polygon( |
| 116 | + data = world, aes(x = x, y = y, group = group), |
| 117 | + fill = LAND_BG, color = NA |
| 118 | + ) + |
| 119 | + # Country borders / coastlines |
| 120 | + geom_path( |
| 121 | + data = world, aes(x = x, y = y, group = group), |
| 122 | + color = INK_SOFT, linewidth = 0.10, alpha = 0.55 |
| 123 | + ) + |
| 124 | + # Standard graticule (30° grid) |
| 125 | + geom_path( |
| 126 | + data = graticule, aes(x = x, y = y, group = group), |
| 127 | + color = INK_SOFT, linewidth = 0.18, alpha = 0.45 |
| 128 | + ) + |
| 129 | + # Equator and prime meridian — slightly bolder |
| 130 | + geom_path( |
| 131 | + data = equator, aes(x = x, y = y, group = group), |
| 132 | + color = INK_SOFT, linewidth = 0.35, alpha = 0.70 |
| 133 | + ) + |
| 134 | + geom_path( |
| 135 | + data = prime_merid, aes(x = x, y = y, group = group), |
| 136 | + color = INK_SOFT, linewidth = 0.35, alpha = 0.70 |
| 137 | + ) + |
| 138 | + # Tissot indicatrices — anyplot brand green (position 1) |
| 139 | + geom_polygon( |
| 140 | + data = tissot, aes(x = x, y = y, group = group), |
| 141 | + fill = ANYPLOT_PALETTE[1], color = PAGE_BG, |
| 142 | + linewidth = 0.06, alpha = 0.78 |
| 143 | + ) + |
| 144 | + # Ellipse outline |
| 145 | + geom_path( |
| 146 | + data = boundary, aes(x = x, y = y), |
| 147 | + color = INK_SOFT, linewidth = 0.40 |
| 148 | + ) + |
| 149 | + coord_fixed() + |
| 150 | + labs( |
| 151 | + title = "map-projections · r · ggplot2 · anyplot.ai", |
| 152 | + subtitle = "Mollweide equal-area projection · Tissot indicatrices (green) confirm equal-area: every circle covers identical surface" |
| 153 | + ) + |
| 154 | + theme_void(base_size = 8) + |
| 155 | + theme( |
| 156 | + plot.background = element_rect(fill = PAGE_BG, color = PAGE_BG), |
| 157 | + plot.title = element_text( |
| 158 | + color = INK, size = 12, hjust = 0.5, |
| 159 | + margin = margin(t = 14, b = 6) |
| 160 | + ), |
| 161 | + plot.subtitle = element_text( |
| 162 | + color = INK_SOFT, size = 9.5, hjust = 0.5, |
| 163 | + margin = margin(b = 12) |
| 164 | + ), |
| 165 | + plot.margin = margin(8, 50, 12, 50) |
| 166 | + ) |
| 167 | + |
| 168 | +# --- Save -------------------------------------------------------------------- |
| 169 | +ggsave( |
| 170 | + filename = sprintf("plot-%s.png", THEME), |
| 171 | + plot = p, |
| 172 | + device = ragg::agg_png, |
| 173 | + width = 8, |
| 174 | + height = 4.5, |
| 175 | + units = "in", |
| 176 | + dpi = 400 |
| 177 | +) |
0 commit comments