Render in custom MapLibre projections (addProjection)

README.md

@@ -4,7 +4,7 @@
 [![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](https://opensource.org/licenses/MIT)
 
 Custom layer for rendering Zarr datasets in MapLibre or Mapbox GL, inspired (and borrowing significant code and concepts from)
-[zarr-gl](https://github.com/carderne/zarr-gl), [zarr-cesium](https://github.com/NOC-OI/zarr-cesium), [@carbonplan/maps](https://github.com/carbonplan/maps), and [deck-gl-raster](https://github.com/developmentseed/deck.gl-raster). Uses [CustomLayerInterface](https://maplibre.org/maplibre-gl-js/docs/API/interfaces/CustomLayerInterface/) to render data directly to the map and supports rendering to globe and mercator projections for both MapLibre and Mapbox. Mapbox's other projections (Albers, Equal Earth, etc.) don't render custom layers, so the layer isn't drawn in them. Input data are reprojected on the fly.
+[zarr-gl](https://github.com/carderne/zarr-gl), [zarr-cesium](https://github.com/NOC-OI/zarr-cesium), [@carbonplan/maps](https://github.com/carbonplan/maps), and [deck-gl-raster](https://github.com/developmentseed/deck.gl-raster). Uses [CustomLayerInterface](https://maplibre.org/maplibre-gl-js/docs/API/interfaces/CustomLayerInterface/) to render data directly to the map and supports rendering to globe and mercator projections for both MapLibre and Mapbox. Mapbox's other projections (Albers, Equal Earth, etc.) don't render custom layers, so the layer isn't drawn in them. On MapLibre the layer also renders in custom map projections registered with `addProjection`, experimentally (see [custom map projections](#custom-map-projections)). Input data are reprojected on the fly.
 
 ## demo
 
@@ -65,6 +65,30 @@ The layer drapes over 3D terrain on both providers, with no configuration needed
 
 With terrain enabled the layer always renders through the draped Mercator-tile path, on the globe too, so the poles show the usual Web Mercator holes on both providers. Each draped tile reads the pyramid level matching its own zoom, so a pitched view over terrain gets fine data near the camera and coarse data at the horizon rather than nothing.
 
+## custom map projections
+
+On MapLibre (maplibre-gl 6.15 or later), the layer renders in planar map projections registered with MapLibre's [`addProjection`](https://maplibre.org/maplibre-gl-js/docs/API/functions/addProjection/), such as Antarctic Polar Stereographic (EPSG:3031). Both that API and this support are experimental and may change in a minor release. Pass the layer the same projection objects given to `addProjection`: it places the data in the projection's world, selects the regions each of the map's tiles covers, and reads the pyramid level matching the projection's own scale.
+
+```ts
+import proj4 from 'proj4'
+import { addProjection, type CRS } from 'maplibre-gl'
+
+const antarctic: CRS = {
+  name: 'EPSG:3031',
+  projection: proj4(
+    'EPSG:4326',
+    '+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
+  ),
+  tileMatrixSet: { origin: [-4194304, 4194304], extentAtZoom0: 8388608 },
+}
+addProjection(antarctic)
+map.setProjection({ type: 'EPSG:3031' })
+
+const layer = new ZarrLayer({ ...options, mapProjections: [antarctic] })
+```
+
+Data in any CRS the layer reads is reprojected into the map's projection. A map in a registered projection draws a single world without copies, and its basemap and terrain tiles must be cut in the same projection. The demo has polar stereographic (EPSG:3031, EPSG:3413), plate carrée, Equal Earth and CONUS Albers (EPSG:5070) maps over Natural Earth GeoJSON, which MapLibre projects in its workers. While the map is in a registered projection that is not among `mapProjections`, the layer is hidden and logs a warning.
+
 ## install
 
 ```bash
@@ -135,6 +159,7 @@ map.on('load', () => {
 | transformRequest | function | - | Transform request URLs and add headers/credentials (see [authentication](#authentication)) |
 | onAuthError | function | - | Called with the HTTP status when a signed request fails with expired credentials (see [authentication](#authentication)) |
 | renderPoles | boolean | `false` | Enable polar coverage in Mapbox globe for EPSG:4326/proj4 datasets (see [globe rendering](#globe-rendering-and-polar-coverage)). No effect on EPSG:3857 data. MapLibre always renders to the poles. |
+| mapProjections | MapProjection[] | `[]` | Projections registered with MapLibre's experimental `addProjection` that the layer renders in: the same objects given to `addProjection` (see [custom map projections](#custom-map-projections)). |
 
 ## methods
 

demo/components/controls.tsx

@@ -14,7 +14,8 @@ import { useThemedColormap } from '@carbonplan/colormaps'
 import { Info, RotatingArrow } from '@carbonplan/icons'
 import { Box, Flex, IconButton } from 'theme-ui'
 import { SidebarDivider } from '@carbonplan/layouts'
-import { useAppStore } from '../lib/store'
+import { useAppStore, type DemoProjection } from '../lib/store'
+import { CUSTOM_PROJECTION_IDS } from '../lib/projections'
 import { subheadingSx } from './shared-controls'
 import DatasetBrowser from './dataset-browser'
 import type {
@@ -306,7 +307,7 @@ const Controls = () => {
   const opacity = useAppStore((state) => state.opacity)
   const clim = useAppStore((state) => state.clim)
   const colormap = useAppStore((state) => state.colormap)
-  const globeProjection = useAppStore((state) => state.globeProjection)
+  const projection = useAppStore((state) => state.projection)
   const terrainEnabled = useAppStore((state) => state.terrainEnabled)
   const renderPoles = useAppStore((state) => state.renderPoles)
   const setRenderPoles = useAppStore((state) => state.setRenderPoles)
@@ -352,7 +353,7 @@ const Controls = () => {
   const setOpacity = useAppStore((state) => state.setOpacity)
   const setClim = useAppStore((state) => state.setClim)
   const setColormap = useAppStore((state) => state.setColormap)
-  const setGlobeProjection = useAppStore((state) => state.setGlobeProjection)
+  const setProjection = useAppStore((state) => state.setProjection)
   const setTerrainEnabled = useAppStore((state) => state.setTerrainEnabled)
   const setMapProvider = useAppStore((state) => state.setMapProvider)
   const setRegionResult = useAppStore((state) => state.setRegionResult)
@@ -870,10 +871,17 @@ const Controls = () => {
         </Column>
         <Column start={2} width={3}>
           <Filter
-            values={{ globe: globeProjection, mercator: !globeProjection }}
+            // The options differ by provider, and Filter expects a fixed set.
+            key={mapProvider}
+            values={Object.fromEntries(
+              (mapProvider === 'maplibre'
+                ? ['globe', 'mercator', ...CUSTOM_PROJECTION_IDS]
+                : ['globe', 'mercator']
+              ).map((option) => [option, projection === option])
+            )}
             setValues={(obj: Record<string, boolean>) => {
-              if (obj.mercator) setGlobeProjection(false)
-              if (obj.globe) setGlobeProjection(true)
+              const selected = Object.keys(obj).find((option) => obj[option])
+              if (selected) setProjection(selected as DemoProjection)
             }}
           />
         </Column>

demo/components/custom-projections.ts

@@ -0,0 +1,234 @@
+import proj4 from 'proj4'
+import {
+  addProjection,
+  importScriptInWorkers,
+  type CRS,
+  type Map as MapLibreMap,
+} from 'maplibre-gl'
+import { isCustomProjection, type CustomProjectionId } from '../lib/projections'
+
+/**
+ * Map projections registered with MapLibre's experimental `addProjection`
+ * (maplibre-gl >= 6.15), over Natural Earth land and borders drawn as GeoJSON.
+ */
+interface CustomProjection {
+  crs: CRS
+  view: { center: [number, number]; zoom: number }
+}
+
+const POLAR_HALF_EXTENT = 4194304
+
+const polar = (
+  name: string,
+  definition: string,
+  center: [number, number]
+): CustomProjection => ({
+  crs: {
+    name,
+    projection: proj4('EPSG:4326', definition),
+    tileMatrixSet: {
+      origin: [-POLAR_HALF_EXTENT, POLAR_HALF_EXTENT],
+      extentAtZoom0: 2 * POLAR_HALF_EXTENT,
+    },
+    definition,
+  },
+  view: { center, zoom: 1.2 },
+})
+
+const LONGLAT = '+proj=longlat +datum=WGS84 +no_defs'
+const EQUAL_EARTH = '+proj=eqearth +lon_0=0 +datum=WGS84 +units=m +no_defs'
+const EQUAL_EARTH_HALF_WIDTH = 17243960
+const CONUS_ALBERS =
+  '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'
+
+export const CUSTOM_PROJECTIONS: Record<CustomProjectionId, CustomProjection> =
+  {
+    antarctic: polar(
+      'EPSG:3031',
+      '+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
+      [0, -90]
+    ),
+    arctic: polar(
+      'EPSG:3413',
+      '+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
+      [0, 90]
+    ),
+    'plate carrée': {
+      crs: {
+        name: 'EPSG:4326',
+        projection: proj4('EPSG:4326', LONGLAT),
+        tileMatrixSet: { origin: [-180, 180], extentAtZoom0: 360 },
+        definition: LONGLAT,
+      },
+      view: { center: [0, 0], zoom: 0.5 },
+    },
+    'equal earth': {
+      crs: {
+        name: 'ESRI:54035',
+        projection: proj4('EPSG:4326', EQUAL_EARTH),
+        tileMatrixSet: {
+          origin: [-EQUAL_EARTH_HALF_WIDTH, EQUAL_EARTH_HALF_WIDTH],
+          extentAtZoom0: 2 * EQUAL_EARTH_HALF_WIDTH,
+        },
+        definition: EQUAL_EARTH,
+      },
+      view: { center: [0, 0], zoom: 0.7 },
+    },
+    'conus albers': {
+      crs: {
+        name: 'EPSG:5070',
+        projection: proj4('EPSG:4326', CONUS_ALBERS),
+        tileMatrixSet: { origin: [-3000000, 4750000], extentAtZoom0: 6000000 },
+        definition: CONUS_ALBERS,
+      },
+      view: { center: [-97, 38], zoom: 1 },
+    },
+  }
+
+/** Passed to the layer's `mapProjections`, the same objects `addProjection` gets. */
+export const MAP_PROJECTIONS: CRS[] = Object.values(CUSTOM_PROJECTIONS).map(
+  (p) => p.crs
+)
+
+const BASEMAP_PREFIX = 'custom-basemap'
+const NATURAL_EARTH =
+  'https://cdn.jsdelivr.net/gh/nvkelso/natural-earth-vector@v5.1.2/geojson'
+const LAND_COLOR = '#2a2e35'
+const BORDER_COLOR = '#5b6069'
+
+let registered = false
+let projectionHandler: Promise<void> | null = null
+/** Bumped on each projection change, so a basemap still loading can tell it is stale. */
+const basemapGenerations = new WeakMap<MapLibreMap, number>()
+
+function register() {
+  if (registered) return
+  registered = true
+  for (const { crs } of Object.values(CUSTOM_PROJECTIONS)) addProjection(crs)
+}
+
+/** GeoJSON sources project in the workers, which build proj4 from `definition`. */
+function loadProjectionHandler(): Promise<void> {
+  projectionHandler ??= importScriptInWorkers(
+    new URL('/projection-handler.mjs', window.location.href).href
+  ).then(
+    () => undefined,
+    (error) => {
+      // Let the next projection change try again.
+      projectionHandler = null
+      throw error
+    }
+  )
+  return projectionHandler
+}
+
+function removeBasemap(map: MapLibreMap) {
+  const style = map.getStyle()
+  for (const layer of style.layers) {
+    if (layer.id.startsWith(BASEMAP_PREFIX)) map.removeLayer(layer.id)
+  }
+  for (const id of Object.keys(style.sources)) {
+    if (id.startsWith(BASEMAP_PREFIX)) map.removeSource(id)
+  }
+}
+
+/** The layer under everything but the background, and the first label. */
+function anchors(map: MapLibreMap) {
+  const { layers } = map.getStyle()
+  return {
+    bottom: layers.find((l) => l.type !== 'background')?.id,
+    labels: layers.find((l) => l.type === 'symbol')?.id,
+  }
+}
+
+function addBasemap(map: MapLibreMap, generation: number) {
+  loadProjectionHandler().then(
+    () => {
+      // The map may have moved on, or gone, while the workers loaded it.
+      if (basemapGenerations.get(map) !== generation) return
+      // Borders draw over the data, under the labels.
+      const { bottom, labels } = anchors(map)
+      const attribution =
+        '<a href="https://www.naturalearthdata.com">Natural Earth</a>'
+      map.addSource(`${BASEMAP_PREFIX}-land`, {
+        type: 'geojson',
+        data: `${NATURAL_EARTH}/ne_50m_land.geojson`,
+        attribution,
+      })
+      // The land polygons are cut at 180°, which would show at the poles.
+      map.addSource(`${BASEMAP_PREFIX}-coast`, {
+        type: 'geojson',
+        data: `${NATURAL_EARTH}/ne_50m_coastline.geojson`,
+      })
+      map.addSource(`${BASEMAP_PREFIX}-borders`, {
+        type: 'geojson',
+        data: `${NATURAL_EARTH}/ne_50m_admin_0_boundary_lines_land.geojson`,
+      })
+      map.addLayer(
+        {
+          id: `${BASEMAP_PREFIX}-land`,
+          type: 'fill',
+          source: `${BASEMAP_PREFIX}-land`,
+          paint: { 'fill-color': LAND_COLOR },
+        },
+        bottom
+      )
+      map.addLayer(
+        {
+          id: `${BASEMAP_PREFIX}-coast`,
+          type: 'line',
+          source: `${BASEMAP_PREFIX}-coast`,
+          paint: { 'line-color': BORDER_COLOR, 'line-width': 0.75 },
+        },
+        labels
+      )
+      map.addLayer(
+        {
+          id: `${BASEMAP_PREFIX}-borders`,
+          type: 'line',
+          source: `${BASEMAP_PREFIX}-borders`,
+          paint: { 'line-color': BORDER_COLOR, 'line-width': 0.5 },
+        },
+        labels
+      )
+    },
+    (error) =>
+      console.warn(
+        '[demo] The GeoJSON basemap could not load its projection handler:',
+        error
+      )
+  )
+}
+
+/**
+ * Put a MapLibre map in one of these projections, or a built-in one, with
+ * the basemap to match: the Web Mercator vector basemap shows only in the
+ * built-in projections.
+ */
+export function setMapProjection(
+  map: MapLibreMap,
+  projection: 'globe' | 'mercator' | CustomProjectionId
+) {
+  const custom = isCustomProjection(projection)
+    ? CUSTOM_PROJECTIONS[projection]
+    : null
+  if (custom) register()
+  if (!basemapGenerations.has(map)) {
+    map.once('remove', () => basemapGenerations.delete(map))
+  }
+  const generation = (basemapGenerations.get(map) ?? 0) + 1
+  basemapGenerations.set(map, generation)
+  const previous = map.getProjection()?.type
+  const next = custom ? custom.crs.name : projection
+  removeBasemap(map)
+  for (const layer of map.getStyle().layers) {
+    if ('source' in layer && layer.source === 'protomaps') {
+      map.setLayoutProperty(layer.id, 'visibility', custom ? 'none' : 'visible')
+    }
+  }
+  map.setProjection({ type: next })
+  if (custom) addBasemap(map, generation)
+  const wasCustom = MAP_PROJECTIONS.some((crs) => crs.name === previous)
+  if (custom && previous !== next) map.jumpTo(custom.view)
+  else if (!custom && wasCustom) map.jumpTo({ center: [0, 20], zoom: 2 })
+}