diff --git a/README.md b/README.md index 80bd663a..331032a4 100644 --- a/README.md +++ b/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 diff --git a/demo/components/controls.tsx b/demo/components/controls.tsx index 6d483fe0..e60e6464 100644 --- a/demo/components/controls.tsx +++ b/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 = () => { [option, projection === option]) + )} setValues={(obj: Record) => { - 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) }} /> diff --git a/demo/components/custom-projections.ts b/demo/components/custom-projections.ts new file mode 100644 index 00000000..1137fd3b --- /dev/null +++ b/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 = + { + 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 | null = null +/** Bumped on each projection change, so a basemap still loading can tell it is stale. */ +const basemapGenerations = new WeakMap() + +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 { + 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 = + 'Natural Earth' + 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 }) +} diff --git a/demo/components/map-shared.tsx b/demo/components/map-shared.tsx index cfb082ee..09620754 100644 --- a/demo/components/map-shared.tsx +++ b/demo/components/map-shared.tsx @@ -10,7 +10,9 @@ import { Map as MapLibreMap, addProtocol, setWorkerUrl } from 'maplibre-gl' import mapboxgl from 'mapbox-gl' import { layers, namedFlavor } from '@protomaps/basemaps' import { Protocol } from 'pmtiles' -import { useAppStore } from '../lib/store' +import { useAppStore, type DemoProjection } from '../lib/store' +import { MAP_PROJECTIONS, setMapProjection } from './custom-projections' +import { isCustomProjection } from '../lib/projections' import type { LayerProps } from '../datasets/types' import MapZoomControls, { useAttributionStyles } from './map-controls' @@ -82,16 +84,16 @@ const mapLibreTheme = { export interface MapConfig { createMap: ( container: HTMLDivElement, - globeProjection: boolean + projection: DemoProjection ) => MapInstance - setProjection: (map: MapInstance, globeProjection: boolean) => void + setProjection: (map: MapInstance, projection: DemoProjection) => void getLayerBeforeId: (map: MapInstance) => string | undefined /** Id of the raster-dem source the terrain toggle drapes onto. */ terrainSource: string } const mapLibreConfig: MapConfig = { - createMap: (container: HTMLDivElement, globeProjection: boolean) => { + createMap: (container: HTMLDivElement, projection: DemoProjection) => { setWorkerUrl( new URL( 'maplibre-gl/dist/maplibre-gl-worker.mjs', @@ -104,7 +106,8 @@ const mapLibreConfig: MapConfig = { return new MapLibreMap({ container, style: { - projection: globeProjection ? { type: 'globe' } : { type: 'mercator' }, + projection: + projection === 'mercator' ? { type: 'mercator' } : { type: 'globe' }, version: 8, glyphs: 'https://carbonplan-maps.s3.us-west-2.amazonaws.com/basemaps/fonts/{fontstack}/{range}.pbf', @@ -133,17 +136,15 @@ const mapLibreConfig: MapConfig = { zoom: window.innerWidth < 640 ? 1.2 : 2.4, }) as MapInstance }, - setProjection: (map: MapInstance, globeProjection: boolean) => { - ;(map as MapLibreMap).setProjection( - globeProjection ? { type: 'globe' } : { type: 'mercator' } - ) + setProjection: (map: MapInstance, projection: DemoProjection) => { + setMapProjection(map as MapLibreMap, projection) }, getLayerBeforeId: () => 'landuse_pedestrian', terrainSource: 'terrain-dem', } const mapboxConfig: MapConfig = { - createMap: (container: HTMLDivElement, globeProjection: boolean) => { + createMap: (container: HTMLDivElement, projection: DemoProjection) => { if (process.env.NEXT_PUBLIC_MAPBOX_ACCESS_TOKEN) { mapboxgl.accessToken = process.env.NEXT_PUBLIC_MAPBOX_ACCESS_TOKEN } @@ -153,7 +154,7 @@ const mapboxConfig: MapConfig = { style: 'mapbox://styles/mapbox/dark-v11', center: [0, 20], zoom: window.innerWidth < 640 ? 1 : 2, - projection: globeProjection ? 'globe' : 'mercator', + projection: projection === 'mercator' ? 'mercator' : 'globe', }) map.on('load', () => { @@ -167,9 +168,9 @@ const mapboxConfig: MapConfig = { return map as MapInstance }, - setProjection: (map: MapInstance, globeProjection: boolean) => { + setProjection: (map: MapInstance, projection: DemoProjection) => { ;(map as mapboxgl.Map).setProjection({ - name: globeProjection ? 'globe' : 'mercator', + name: projection === 'mercator' ? 'mercator' : 'globe', } as mapboxgl.ProjectionSpecification) }, getLayerBeforeId: (map: MapInstance) => { @@ -241,6 +242,7 @@ export const useMapLayer = (map: MapInstance | null, isMapLoaded: boolean) => { proj4: datasetModule.proj4, onLoadingStateChange: setLoadingState, renderPoles, + mapProjections: MAP_PROJECTIONS, } if (datasetModule.store) { @@ -411,7 +413,7 @@ export const Map = () => { const sidebarWidth = useAppStore((state) => state.sidebarWidth) const mapProvider = useAppStore((state) => state.mapProvider) - const globeProjection = useAppStore((state) => state.globeProjection) + const projection = useAppStore((state) => state.projection) const terrainEnabled = useAppStore((state) => state.terrainEnabled) const loadingState = useAppStore((state) => state.loadingState) const setMapInstance = useAppStore((state) => state.setMapInstance) @@ -421,7 +423,7 @@ export const Map = () => { useEffect(() => { if (!mapContainer.current) return - const newMap = mapConfig.createMap(mapContainer.current, globeProjection) + const newMap = mapConfig.createMap(mapContainer.current, projection) mapInstanceRef.current = newMap newMap.on('load', () => { @@ -445,13 +447,15 @@ export const Map = () => { useEffect(() => { if (!map || !isMapLoaded) return - mapConfig.setProjection(map, globeProjection) - }, [map, isMapLoaded, globeProjection]) + mapConfig.setProjection(map, projection) + }, [map, isMapLoaded, projection]) useEffect(() => { if (!map || !isMapLoaded) return try { - if (terrainEnabled) { + // The terrain tiles are Web Mercator, so they cannot drape the grid of + // a custom projection. + if (terrainEnabled && !isCustomProjection(projection)) { map.setTerrain({ source: mapConfig.terrainSource, exaggeration: 1.5 }) } else { map.setTerrain(null) @@ -459,7 +463,7 @@ export const Map = () => { } catch (e) { console.warn('Error toggling terrain:', e) } - }, [map, isMapLoaded, terrainEnabled, mapConfig]) + }, [map, isMapLoaded, terrainEnabled, projection, mapConfig]) useMapLayer(map, isMapLoaded) diff --git a/demo/lib/projections.ts b/demo/lib/projections.ts new file mode 100644 index 00000000..f04f5ca3 --- /dev/null +++ b/demo/lib/projections.ts @@ -0,0 +1,15 @@ +/** Map projections the demo registers with MapLibre's `addProjection`. */ +export const CUSTOM_PROJECTION_IDS = [ + 'antarctic', + 'arctic', + 'plate carrée', + 'equal earth', + 'conus albers', +] as const +export type CustomProjectionId = (typeof CUSTOM_PROJECTION_IDS)[number] + +export function isCustomProjection( + projection: string +): projection is CustomProjectionId { + return (CUSTOM_PROJECTION_IDS as readonly string[]).includes(projection) +} diff --git a/demo/lib/store.ts b/demo/lib/store.ts index 5745ee2d..1c174a55 100644 --- a/demo/lib/store.ts +++ b/demo/lib/store.ts @@ -1,6 +1,7 @@ import { create } from 'zustand' import { DATASET_MAP, DEFAULT_DATASET_ID } from '../datasets' import type { Dataset } from '../datasets/types' +import { isCustomProjection, type CustomProjectionId } from './projections' import type { MapProvider, MapInstance } from '../components/map-shared' import type { ZarrLayer, @@ -8,13 +9,15 @@ import type { LoadingState, } from '@carbonplan/zarr-layer' +export type DemoProjection = 'globe' | 'mercator' | CustomProjectionId + interface AppState { sidebarWidth: number | null datasetId: string opacity: number clim: [number, number] colormap: string - globeProjection: boolean + projection: DemoProjection terrainEnabled: boolean renderPoles: boolean mapProvider: MapProvider @@ -30,7 +33,7 @@ interface AppState { setOpacity: (opacity: number) => void setClim: (clim: [number, number]) => void setColormap: (colormap: string) => void - setGlobeProjection: (globeProjection: boolean) => void + setProjection: (projection: DemoProjection) => void setTerrainEnabled: (terrainEnabled: boolean) => void setRenderPoles: (renderPoles: boolean) => void setMapProvider: (provider: MapProvider) => void @@ -52,7 +55,7 @@ export const useAppStore = create((set, get) => ({ opacity: 1, clim: defaultModule.clim, colormap: defaultModule.colormap, - globeProjection: true, + projection: 'globe', terrainEnabled: false, renderPoles: false, mapProvider: 'maplibre', @@ -79,10 +82,18 @@ export const useAppStore = create((set, get) => ({ setOpacity: (opacity) => set({ opacity }), setClim: (clim) => set({ clim }), setColormap: (colormap) => set({ colormap }), - setGlobeProjection: (globeProjection) => set({ globeProjection }), + setProjection: (projection) => set({ projection }), setTerrainEnabled: (terrainEnabled) => set({ terrainEnabled }), setRenderPoles: (renderPoles) => set({ renderPoles }), - setMapProvider: (mapProvider) => set({ mapProvider }), + // Custom projections are MapLibre's addProjection, which Mapbox lacks. + setMapProvider: (mapProvider) => + set((state) => ({ + mapProvider, + projection: + mapProvider === 'mapbox' && isCustomProjection(state.projection) + ? 'globe' + : state.projection, + })), setActiveDatasetState: (updates) => { set((state) => ({ datasetState: { ...state.datasetState, ...updates }, diff --git a/demo/package-lock.json b/demo/package-lock.json index 29d0055c..c2a98b84 100644 --- a/demo/package-lock.json +++ b/demo/package-lock.json @@ -17,9 +17,10 @@ "@protomaps/basemaps": "^5.0.1", "icechunk-js": "^0.6.1", "mapbox-gl": "3.32.0", - "maplibre-gl": "^6.14.0", + "maplibre-gl": "^6.16.0", "next": "^16.4.2", "pmtiles": "^4.5.0", + "proj4": "^2.22.0", "react": "^19.3.1", "react-dom": "^19.3.1", "theme-ui": "^0.17.4", @@ -1085,9 +1086,9 @@ } }, "node_modules/@maplibre/geojson-vt": { - "version": "6.2.0", - "resolved": "https://registry.npmjs.org/@maplibre/geojson-vt/-/geojson-vt-6.2.0.tgz", - "integrity": "sha512-r24Cj+32JFQG4tscFoXpVreNJRqMmP2PxDNRaU89SyDJKKC3SWTEuEUgf76JhGdR7uRv6LAr7bcOE4GWWUfmGw==", + "version": "6.3.0", + "resolved": "https://registry.npmjs.org/@maplibre/geojson-vt/-/geojson-vt-6.3.0.tgz", + "integrity": "sha512-3Zw2jsU9pIJG3G2CQFbLjWj2MlAHodhD7cSY+15iywUuLmmkgzcKRnueVRZ4AcC/5wOtfhidcJAAzBpoBNs3Xg==", "license": "ISC", "dependencies": { "kdbush": "^4.1.0" @@ -1122,14 +1123,14 @@ } }, "node_modules/@maplibre/vt-pbf": { - "version": "4.3.2", - "resolved": "https://registry.npmjs.org/@maplibre/vt-pbf/-/vt-pbf-4.3.2.tgz", - "integrity": "sha512-j6p0AdjvAR19Z3XaCysle7A4ZSo08tYOzxD0Y9NQylwPAkwJJeYub5b2eVucdeDh7erhv69DahoLOevDRERRUw==", + "version": "4.3.3", + "resolved": "https://registry.npmjs.org/@maplibre/vt-pbf/-/vt-pbf-4.3.3.tgz", + "integrity": "sha512-dEe5XUd64RQbxNxy2OzIkJnwNDWjVmSXjCd0A7zLirm7QJ/WWbmtDFIkfsiM4yt1cplqzG6oSivRYeWx1SM1yg==", "license": "MIT", "dependencies": { "@mapbox/point-geometry": "^1.1.0", "@types/geojson": "^7946.0.16", - "pbf": "^5.1.0" + "pbf": "^5.1.2" } }, "node_modules/@msgpack/msgpack": { @@ -1992,19 +1993,19 @@ ] }, "node_modules/maplibre-gl": { - "version": "6.14.0", - "resolved": "https://registry.npmjs.org/maplibre-gl/-/maplibre-gl-6.14.0.tgz", - "integrity": "sha512-L7U9y9TyPJ3P/JJbW+eH+A3FzdPq+3fnZFqqeCwf/S6cRGganHwVI5b/C8/zOYzreUzsloDPZVVaMQA94wOt4A==", + "version": "6.16.0", + "resolved": "https://registry.npmjs.org/maplibre-gl/-/maplibre-gl-6.16.0.tgz", + "integrity": "sha512-gApO66EMtoDZuio680RUBvskR7ahFd/NNhLN4tBAs+Y7FeT4gQV8zVFgYVMRmFxQznOEzecW19w6xOd9XjZjiA==", "license": "BSD-3-Clause", "dependencies": { "@mapbox/point-geometry": "^1.1.0", - "@mapbox/tiny-sdf": "^2.2.0", + "@mapbox/tiny-sdf": "^2.3.0", "@mapbox/unitbezier": "^1.0.0", "@mapbox/vector-tile": "^3.0.0", - "@maplibre/geojson-vt": "^6.1.2", + "@maplibre/geojson-vt": "^6.3.0", "@maplibre/maplibre-gl-style-spec": "^26.4.4", "@maplibre/mlt": "^1.3.0", - "@maplibre/vt-pbf": "^4.3.2", + "@maplibre/vt-pbf": "^4.3.3", "@types/geojson": "^7946.0.16", "bidi-js": "^1.1.0", "earcut": "^3.2.4", @@ -2024,6 +2025,12 @@ "url": "https://github.com/maplibre/maplibre-gl-js?sponsor=1" } }, + "node_modules/mgrs": { + "version": "1.0.0", + "resolved": "https://registry.npmjs.org/mgrs/-/mgrs-1.0.0.tgz", + "integrity": "sha512-awNbTOqCxK1DBGjalK3xqWIstBZgN6fxsMSiXLs9/spqWkF2pAhb2rrYCFSsr1/tT7PhcDGjZndG8SWYn0byYA==", + "license": "MIT" + }, "node_modules/minimist": { "version": "1.2.8", "resolved": "https://registry.npmjs.org/minimist/-/minimist-1.2.8.tgz", @@ -2286,6 +2293,27 @@ "url": "https://github.com/prettier/prettier?sponsor=1" } }, + "node_modules/proj4": { + "version": "2.22.0", + "resolved": "https://registry.npmjs.org/proj4/-/proj4-2.22.0.tgz", + "integrity": "sha512-MVvH09QrYm3A235rcX39WDiCKanrBiblu28l0KYIn/B2DofmYhvP/LLS43kE66jNDfx/P1l/nah3JmsAq8UgbQ==", + "license": "MIT", + "dependencies": { + "mgrs": "1.0.0", + "wkt-parser": "^1.5.5" + }, + "funding": { + "url": "https://github.com/sponsors/ahocevar" + }, + "peerDependencies": { + "geotiff": "*" + }, + "peerDependenciesMeta": { + "geotiff": { + "optional": true + } + } + }, "node_modules/protocol-buffers-schema": { "version": "3.6.1", "resolved": "https://registry.npmjs.org/protocol-buffers-schema/-/protocol-buffers-schema-3.6.1.tgz", @@ -2613,6 +2641,15 @@ "node": ">=18" } }, + "node_modules/wkt-parser": { + "version": "1.6.0", + "resolved": "https://registry.npmjs.org/wkt-parser/-/wkt-parser-1.6.0.tgz", + "integrity": "sha512-O2WMUagkdOKkgcfqZBUauysHVeifCwcs+zUd2QG7UIzVljLwK9g5xHUfL/TkQFwmQ3eF7oHFKGnMcb/kZdvbaw==", + "license": "MIT", + "funding": { + "url": "https://github.com/sponsors/ahocevar" + } + }, "node_modules/yaml": { "version": "1.10.3", "resolved": "https://registry.npmjs.org/yaml/-/yaml-1.10.3.tgz", diff --git a/demo/package.json b/demo/package.json index 05a1e5a9..4177fef1 100644 --- a/demo/package.json +++ b/demo/package.json @@ -22,9 +22,10 @@ "@protomaps/basemaps": "^5.0.1", "icechunk-js": "^0.6.1", "mapbox-gl": "3.32.0", - "maplibre-gl": "^6.14.0", + "maplibre-gl": "^6.16.0", "next": "^16.4.2", "pmtiles": "^4.5.0", + "proj4": "^2.22.0", "react": "^19.3.1", "react-dom": "^19.3.1", "theme-ui": "^0.17.4", diff --git a/demo/public/projection-handler.mjs b/demo/public/projection-handler.mjs new file mode 100644 index 00000000..9e0911a4 --- /dev/null +++ b/demo/public/projection-handler.mjs @@ -0,0 +1,5 @@ +// Loaded into MapLibre's workers so GeoJSON sources follow a registered +// projection: the workers build its converter from the CRS's `definition`. +import proj4 from 'https://cdn.jsdelivr.net/npm/proj4@2.22.0/+esm' + +self.setProjectionHandler((definition) => proj4(definition)) diff --git a/src/clip-boxes.test.ts b/src/clip-boxes.test.ts index 2820fc7e..e9a03916 100644 --- a/src/clip-boxes.test.ts +++ b/src/clip-boxes.test.ts @@ -58,7 +58,8 @@ describe('tileClipPasses', () => { const rows = absolute(pass, mesh) expect(rows).toHaveLength(3) for (const row of rows) { - expect(row.x0).toBeCloseTo(0, 6) + // The rows start on the mesh's west edge, which opens. + expect(row.x0).toBeLessThan(0) expect(row.x1).toBeCloseTo(0.5, 6) } @@ -68,6 +69,24 @@ describe('tileClipPasses', () => { expect(passes.map((p) => p.length / 4)).toEqual([MAX_CLIP_BOXES, 8]) }) + it('opens box edges that lie on the mesh edge, so neighbours leave no gap', () => { + // A region whose east edge is a tile edge: the box there would decide + // the region's edge pixels by a rounding error against the delta the + // shader interpolates, so it moves out a mesh width. Its west edge, + // inside the mesh, still clips. + const owner = tile(2, 1, 1) + const bounds = getTileBounds(owner.tileId) + const mesh = { x0: 0.2, y0: 0.2, x1: bounds.x1, y1: 0.6 } + + const [pass] = tileClipPasses(mesh, [owner]) + + const [box] = absolute(pass, mesh) + expect(box.x0).toBeCloseTo(bounds.x0, 6) + expect(box.y0).toBeCloseTo(bounds.y0, 6) + expect(box.y1).toBeCloseTo(bounds.y1, 6) + expect(box.x1).toBeCloseTo(mesh.x1 + (mesh.x1 - mesh.x0), 6) + }) + it('meets tiles across the antimeridian through the next world copy', () => { // A region unwrapped past 180 degrees spans x = 0.99..1.01. The tile on // the far side lives at x = 0..1/64 in the base world; its copy one diff --git a/src/clip-boxes.ts b/src/clip-boxes.ts index 9afe83b3..76ee666f 100644 --- a/src/clip-boxes.ts +++ b/src/clip-boxes.ts @@ -49,6 +49,28 @@ export function tileClipBounds({ tileId, poles }: ClipTile): Box { return bounds } +/** Box and mesh edges closer than this are the same edge. */ +const SHARED_EDGE = 1e-12 + +/** + * A box edge on the mesh's own edge clips nothing the mesh draws, but the + * compare there falls either way by a rounding error: where a tile edge + * meets a region edge, the region's edge pixels can fail it while the + * neighbouring region does not reach them, leaving a hairline gap. Move + * such edges a mesh width out, so only edges inside the mesh clip. + */ +function openMeshEdges(box: Box, mesh: MeshMercatorBounds): Box { + const width = mesh.x1 - mesh.x0 + const height = mesh.y1 - mesh.y0 + const on = (a: number, b: number) => Math.abs(a - b) <= SHARED_EDGE + return { + x0: on(box.x0, mesh.x0) ? mesh.x0 - width : box.x0, + y0: on(box.y0, mesh.y0) ? mesh.y0 - height : box.y0, + x1: on(box.x1, mesh.x1) ? mesh.x1 + width : box.x1, + y1: on(box.y1, mesh.y1) ? mesh.y1 + height : box.y1, + } +} + const overlaps = (a: Box, b: Box) => a.x0 < b.x1 && a.x1 > b.x0 && a.y0 < b.y1 && a.y1 > b.y0 @@ -115,7 +137,7 @@ export function tileClipPasses( if (overlaps(box, mesh)) boxes.push(box) } } - const merged = mergeRows(boxes) + const merged = mergeRows(boxes).map((box) => openMeshEdges(box, mesh)) const passes: Float32Array[] = [] for (let start = 0; start < merged.length; start += MAX_CLIP_BOXES) { const chunk = merged.slice(start, start + MAX_CLIP_BOXES) diff --git a/src/index.ts b/src/index.ts index 01a3d8fe..b1d70ae3 100644 --- a/src/index.ts +++ b/src/index.ts @@ -12,6 +12,7 @@ export type { RequestParameters, } from './types' export type { UniformValue } from './renderer-types' +export type { MapProjection } from './world-mapping' // Query interface exports export type { diff --git a/src/map-utils.projection.test.ts b/src/map-utils.projection.test.ts index cb24b8c9..0b5c6b6b 100644 --- a/src/map-utils.projection.test.ts +++ b/src/map-utils.projection.test.ts @@ -1,6 +1,7 @@ import { describe, it, expect } from 'vitest' import { computeWorldOffsets, + customProjectionName, isGlobeProjection, resolveProjectionParams, } from './map-utils' @@ -54,6 +55,36 @@ describe('isGlobeProjection', () => { }) }) +describe('customProjectionName', () => { + it('names projections registered with addProjection', () => { + expect(customProjectionName({ type: 'EPSG:3031' })).toBe('EPSG:3031') + expect(customProjectionName({ type: 'custom:floor-plan' })).toBe( + 'custom:floor-plan' + ) + expect(customProjectionName({ type: 'identity' })).toBe('identity') + }) + + it('returns null for built-in projections', () => { + expect(customProjectionName({ type: 'mercator' })).toBe(null) + expect(customProjectionName({ type: 'globe' })).toBe(null) + expect(customProjectionName({ type: 'vertical-perspective' })).toBe(null) + expect( + customProjectionName({ + type: [ + 'interpolate', + ['linear'], + ['zoom'], + 10, + 'globe', + 12, + 'mercator', + ], + }) + ).toBe(null) + expect(customProjectionName(null)).toBe(null) + }) +}) + describe('resolveProjectionParams — MapLibre', () => { it('takes the matrix and projection data from the params object', () => { const resolved = resolveProjectionParams(maplibreParams(0)) diff --git a/src/map-utils.ts b/src/map-utils.ts index a68ba98b..8d290bb1 100644 --- a/src/map-utils.ts +++ b/src/map-utils.ts @@ -250,6 +250,19 @@ export function isGlobeProjection( return projection?.type === 'globe' || projection?.name === 'globe' } +/** + * The name of a projection registered with MapLibre's `addProjection` + * (maplibre-gl >= 6.15), or null for a built-in projection. Registered names + * contain a colon; `identity` is the one MapLibre registers itself. + */ +export function customProjectionName( + projection: { type?: unknown } | null | undefined +): string | null { + const type = projection?.type + if (typeof type !== 'string') return null + return type.includes(':') || type === 'identity' ? type : null +} + interface ProjectionResolution { matrix: number[] | Float32Array | Float64Array | null shaderData?: ShaderData diff --git a/src/mesh-reprojector.ts b/src/mesh-reprojector.ts index 14c28d07..e6c03450 100644 --- a/src/mesh-reprojector.ts +++ b/src/mesh-reprojector.ts @@ -855,3 +855,340 @@ export function createHybridMesh( meshBounds, } } + +/** + * How far, as a share of an edge's length, the projected middle of an edge + * may land from the middle of the drawn edge. Past it the edge crosses a seam + * or fold in the projection, such as a world projection's antimeridian, and + * would draw across the map. + */ +const SEAM_MIDPOINT_RATIO = 0.25 +/** + * An edge that strays that far is halved toward its widest part this many + * times, homing in on any jump. A seam leaves a gap the data doesn't fill, so + * the gap's middle maps back far from the edge; a bend, even one sharpened + * without bound like a cone's tip, maps back onto it. + */ +const SEAM_BISECTIONS = 30 +/** Under a tenth of a pixel at zoom 24, in world units. */ +const SEAM_MIN_LENGTH = 1e-11 + +interface WorldMeshOptions { + geoBounds: { xMin: number; xMax: number; yMin: number; yMax: number } + width: number + height: number + lonSubdivisions: number + latSubdivisions: number + /** Source CRS to the map's world square and back. */ + toWorld: (x: number, y: number) => [number, number] + fromWorld: (x: number, y: number) => [number, number] + latIsAscending: boolean + /** Source x is longitude, so positions 360 apart are the same place. */ + allowUnwrappedLongitudes?: boolean + maxError?: number +} + +/** The part of a polygon inside the unit square. */ +function clipToUnitSquare(polygon: [number, number][]): [number, number][] { + let points = polygon + for (const [axis, bound, keepBelow] of [ + [0, 0, false], + [0, 1, true], + [1, 0, false], + [1, 1, true], + ] as const) { + const inside = (p: [number, number]) => + keepBelow ? p[axis] <= bound : p[axis] >= bound + const clipped: [number, number][] = [] + for (let i = 0; i < points.length; i++) { + const p = points[i] + const q = points[(i + 1) % points.length] + if (inside(p)) clipped.push(p) + if (inside(p) !== inside(q)) { + const f = (bound - p[axis]) / (q[axis] - p[axis]) + clipped.push([p[0] + (q[0] - p[0]) * f, p[1] + (q[1] - p[1]) * f]) + } + } + points = clipped + if (points.length === 0) break + } + return points +} + +/** + * Mesh for a map in a projection registered with `addProjection`: vertices + * go straight from the source CRS to the map's world square, refined where + * that mapping bends, and are encoded as deltas from the region's own + * anchor like the Mercator mesh. Triangles with a vertex the map's + * projection cannot place, spanning a seam in it, or wrapping the map + * around a point it sends to infinity are dropped. + */ +export function createWorldMesh(options: WorldMeshOptions): AdaptiveMeshResult { + const { + geoBounds, + width, + height, + lonSubdivisions, + latSubdivisions, + toWorld, + fromWorld, + latIsAscending, + allowUnwrappedLongitudes = false, + maxError = DEFAULT_MESH_MAX_ERROR, + } = options + const { xMin, xMax, yMin, yMax } = geoBounds + + const reprojector = createReprojector({ + bounds: [xMin, yMin, xMax, yMax], + width, + height, + latIsAscending, + transformer: { forward: toWorld, inverse: fromWorld }, + }) + for ( + let i = 0; + i < MAX_ITERATIONS && reprojector.getMaxError() > maxError; + i++ + ) { + const prevVertCount = reprojector.uvs.length / 2 + reprojector.refine() + const newVertCount = reprojector.uvs.length / 2 + if ( + newVertCount >= MAX_ADAPTIVE_VERTICES || + newVertCount === prevVertCount + ) { + break + } + } + + // Edges get a vertex at every data pixel and no others, so neighboring + // regions share the vertices along the edge between them and meet without + // a gap where the projection curves it. + const adaptiveUVs: number[] = [] + for (let k = 0; k < reprojector.uvs.length; k += 2) { + const u = reprojector.uvs[k] + const v = reprojector.uvs[k + 1] + if (u > 0 && u < 1 && v > 0 && v < 1) adaptiveUVs.push(u, v) + } + const cols = Math.max(1, Math.ceil(lonSubdivisions)) + const rows = Math.max(1, Math.ceil(latSubdivisions)) + const edgeCols = Math.max(1, Math.ceil(width)) + const edgeRows = Math.max(1, Math.ceil(height)) + const uvs = new Float64Array( + adaptiveUVs.length + + (cols - 1) * (rows - 1) * 2 + + (edgeCols + 1 + edgeRows + 1) * 4 + ) + uvs.set(adaptiveUVs) + let offset = adaptiveUVs.length + for (let row = 1; row < rows; row++) { + for (let col = 1; col < cols; col++) { + uvs[offset++] = col / cols + uvs[offset++] = row / rows + } + } + for (let col = 0; col <= edgeCols; col++) { + for (const v of [0, 1]) { + uvs[offset++] = col / edgeCols + uvs[offset++] = v + } + } + for (let row = 0; row <= edgeRows; row++) { + for (const u of [0, 1]) { + uvs[offset++] = u + uvs[offset++] = row / edgeRows + } + } + const triangles = new Delaunator(uvs).triangles + + const numVerts = uvs.length / 2 + const world = new Float64Array(numVerts * 2) + let minX = Infinity + let minY = Infinity + let maxX = -Infinity + let maxY = -Infinity + const sourceAt = (u: number, v: number): [number, number] => [ + xMin + u * (xMax - xMin), + latIsAscending ? yMin + v * (yMax - yMin) : yMax - v * (yMax - yMin), + ] + const worldAt = (u: number, v: number) => toWorld(...sourceAt(u, v)) + for (let i = 0; i < numVerts; i++) { + const [x, y] = worldAt(uvs[i * 2], uvs[i * 2 + 1]) + world[i * 2] = x + world[i * 2 + 1] = y + } + + const finite = (i: number) => + Number.isFinite(world[i * 2]) && Number.isFinite(world[i * 2 + 1]) + // Each edge is checked once, though up to two triangles share it. + const seams = new Map() + const hasGap = (a: number, b: number) => { + let u0 = uvs[a * 2] + let v0 = uvs[a * 2 + 1] + let u1 = uvs[b * 2] + let v1 = uvs[b * 2 + 1] + let p0 = [world[a * 2], world[a * 2 + 1]] + let p1 = [world[b * 2], world[b * 2 + 1]] + for (let i = 0; i < SEAM_BISECTIONS; i++) { + const um = (u0 + u1) / 2 + const vm = (v0 + v1) / 2 + const pm = worldAt(um, vm) + if (!Number.isFinite(pm[0]) || !Number.isFinite(pm[1])) return true + const before = Math.hypot(pm[0] - p0[0], pm[1] - p0[1]) + const after = Math.hypot(p1[0] - pm[0], p1[1] - pm[1]) + if (Math.max(before, after) <= SEAM_MIN_LENGTH) return false + if (before >= after) { + u1 = um + v1 = vm + p1 = pm + } else { + u0 = um + v0 = vm + p0 = pm + } + } + const [sx, sy] = fromWorld((p0[0] + p1[0]) / 2, (p0[1] + p1[1]) / 2) + if (!Number.isFinite(sx) || !Number.isFinite(sy)) return true + const [jx, jy] = sourceAt((u0 + u1) / 2, (v0 + v1) / 2) + let dx = sx - jx + if (allowUnwrappedLongitudes) dx -= 360 * Math.round(dx / 360) + return ( + Math.hypot(dx / (xMax - xMin), (sy - jy) / (yMax - yMin)) > + Math.hypot(uvs[a * 2] - uvs[b * 2], uvs[a * 2 + 1] - uvs[b * 2 + 1]) + ) + } + const crossesSeam = (a: number, b: number) => { + const key = Math.min(a, b) * numVerts + Math.max(a, b) + let crosses = seams.get(key) + if (crosses === undefined) { + const [mx, my] = worldAt( + (uvs[a * 2] + uvs[b * 2]) / 2, + (uvs[a * 2 + 1] + uvs[b * 2 + 1]) / 2 + ) + const length = Math.hypot( + world[a * 2] - world[b * 2], + world[a * 2 + 1] - world[b * 2 + 1] + ) + crosses = + !( + Math.hypot( + mx - (world[a * 2] + world[b * 2]) / 2, + my - (world[a * 2 + 1] + world[b * 2 + 1]) / 2 + ) <= + length * SEAM_MIDPOINT_RATIO + ) && hasGap(a, b) + seams.set(key, crosses) + } + return crosses + } + // The map draws nothing past the edge of its tile matrix set. + const outside = (a: number, b: number, c: number) => + [0, 1].some( + (axis) => + Math.max( + world[a * 2 + axis], + world[b * 2 + axis], + world[c * 2 + axis] + ) < 0 || + Math.min( + world[a * 2 + axis], + world[b * 2 + axis], + world[c * 2 + axis] + ) > 1 + ) + // A point the projection sends to infinity, like the far pole of a polar + // stereographic map, can sit inside a triangle whose corners all land off + // the map. The triangle then spans the map though none of its data is + // there, so a triangle with every corner off the map is kept only if a + // point it draws on the map maps back to it. + const offMap = (i: number) => + world[i * 2] < 0 || + world[i * 2] > 1 || + world[i * 2 + 1] < 0 || + world[i * 2 + 1] > 1 + const wrapsMap = (a: number, b: number, c: number) => { + if (!offMap(a) || !offMap(b) || !offMap(c)) return false + const onMap = clipToUnitSquare( + [a, b, c].map((i): [number, number] => [world[i * 2], world[i * 2 + 1]]) + ) + if (onMap.length === 0) return true + let px = 0 + let py = 0 + for (const [x, y] of onMap) { + px += x / onMap.length + py += y / onMap.length + } + let [sx, sy] = fromWorld(px, py) + if (!Number.isFinite(sx) || !Number.isFinite(sy)) return true + const us = [a, b, c].map((i) => uvs[i * 2]) + const vs = [a, b, c].map((i) => uvs[i * 2 + 1]) + if (allowUnwrappedLongitudes) { + const cx = xMin + ((us[0] + us[1] + us[2]) / 3) * (xMax - xMin) + sx += 360 * Math.round((cx - sx) / 360) + } + const u = (sx - xMin) / (xMax - xMin) + const v = latIsAscending + ? (sy - yMin) / (yMax - yMin) + : (yMax - sy) / (yMax - yMin) + // The drawn triangle is straight where the data's image curves, so the + // point may map back beside the triangle, but not further than its size. + const reach = Math.max( + Math.max(...us) - Math.min(...us), + Math.max(...vs) - Math.min(...vs) + ) + return ( + u < Math.min(...us) - reach || + u > Math.max(...us) + reach || + v < Math.min(...vs) - reach || + v > Math.max(...vs) + reach + ) + } + const kept: number[] = [] + for (let t = 0; t < triangles.length; t += 3) { + const a = triangles[t] + const b = triangles[t + 1] + const c = triangles[t + 2] + if (!finite(a) || !finite(b) || !finite(c)) continue + if (outside(a, b, c)) continue + if (crossesSeam(a, b) || crossesSeam(b, c) || crossesSeam(c, a)) continue + if (wrapsMap(a, b, c)) continue + kept.push(a, b, c) + for (const i of [a, b, c]) { + minX = Math.min(minX, world[i * 2]) + minY = Math.min(minY, world[i * 2 + 1]) + maxX = Math.max(maxX, world[i * 2]) + maxY = Math.max(maxY, world[i * 2 + 1]) + } + } + + if (kept.length === 0) { + return { + positions: new Float32Array(0), + texCoords: new Float32Array(0), + indices: new Uint32Array(0), + meshBounds: { x0: 0, y0: 0, x1: 0, y1: 0 }, + } + } + + const anchorX = (minX + maxX) / 2 + const anchorY = (minY + maxY) / 2 + const halfX = (maxX - minX) / 2 || 0.5 + const halfY = (maxY - minY) / 2 || 0.5 + const positions = new Float32Array(numVerts * 2) + for (let i = 0; i < numVerts; i++) { + positions[i * 2] = finite(i) ? (world[i * 2] - anchorX) / halfX : 0 + positions[i * 2 + 1] = finite(i) ? (world[i * 2 + 1] - anchorY) / halfY : 0 + } + + return { + positions, + texCoords: new Float32Array(uvs), + indices: new Uint32Array(kept), + meshBounds: { + x0: anchorX - halfX, + y0: anchorY - halfY, + x1: anchorX + halfX, + y1: anchorY + halfY, + }, + } +} diff --git a/src/region-math.ts b/src/region-math.ts index 84aaf56e..eb452ff3 100644 --- a/src/region-math.ts +++ b/src/region-math.ts @@ -10,6 +10,12 @@ import type { ProjectionContext } from './projection-utils' import { sampleEdgesToMercatorBounds } from './projection-utils' import type { LevelMeta } from './region-state' import type { DimIndicesProps, MapLike, ResolutionLevel } from './types' +import { + worldBoundsOf, + worldToSource, + type WorldBounds, + type WorldMapping, +} from './world-mapping' export type RegionCoordinate = { regionX: number; regionY: number } type SourceBounds = { xMin: number; xMax: number; yMin: number; yMax: number } @@ -342,6 +348,145 @@ export function getCandidateRegions({ return candidates } +/** Below this overlap a region only touches a world box. */ +const STRICT_OVERLAP_WORLD = 1e-9 + +/** + * Regions of a level whose footprint in the world square of a map in a + * registered projection overlaps a box. Candidates come from the box mapped + * back to the source CRS; each is kept only when its own world footprint + * overlaps by more than a touch. + */ +export function getRegionsInWorldBox({ + box, + xyLimits, + levelMeta, + projection, + mapping, + latIsAscending, +}: { + box: WorldBounds + xyLimits: XYLimits | null + levelMeta: LevelMeta | null + projection: ProjectionContext + mapping: WorldMapping + latIsAscending: boolean +}): RegionCoordinate[] { + const to4326 = projection.to4326 + if (!xyLimits || !levelMeta || !to4326) return [] + const { width, height, regionSize } = levelMeta + const [regionH, regionW] = regionSize + const limits = levelMeta.xyLimits ?? xyLimits + const wrapLongitude = + projection.kind === 'epsg4326' && (limits.xMin < -180 || limits.xMax > 180) + const candidates = getCandidateRegions({ + west: box.x0, + south: box.y0, + east: box.x1, + north: box.y1, + transformer: { forward: worldToSource(mapping, to4326) }, + numRegionsX: Math.ceil(width / regionW), + numRegionsY: Math.ceil(height / regionH), + regionW, + regionH, + width, + height, + xyLimits: limits, + latIsAscending, + wrapLongitude, + }) + const inset = STRICT_OVERLAP_WORLD + const reach = sourceReachOfWorldBox(box, mapping, to4326) + // The map gives longitudes in [-180, 180]; regions past it are reached a + // turn away. + const turns = wrapLongitude ? [-360, 0, 360] : [0] + return candidates.filter(({ regionX, regionY }) => { + const bounds = getRegionBounds({ + regionX, + regionY, + levelMeta, + xyLimits, + latIsAscending, + }) + // Curved footprints, such as lon/lat cells around a pole, have boxes far + // larger than themselves; both boxes must overlap. + if ( + reach && + (bounds.yMax <= reach.yMin || + bounds.yMin >= reach.yMax || + turns.every( + (turn) => + bounds.xMax + turn <= reach.xMin || bounds.xMin + turn >= reach.xMax + )) + ) { + return false + } + const footprint = worldBoundsOf(bounds, mapping, to4326, 4) + return ( + !!footprint && + footprint.x0 < box.x1 - inset && + footprint.x1 > box.x0 + inset && + footprint.y0 < box.y1 - inset && + footprint.y1 > box.y0 + inset + ) + }) +} + +/** + * The source-CRS box a world box reaches, from a grid of samples across it + * and each pole inside it, padded by a grid cell for what falls between + * samples. Null when no sample has a source position. + */ +function sourceReachOfWorldBox( + box: WorldBounds, + mapping: WorldMapping, + to4326: { + forward(x: number, y: number): [number, number] + inverse(lng: number, lat: number): [number, number] + } +): SourceBounds | null { + const toSource = worldToSource(mapping, to4326) + const samples = 8 + let xMin = Infinity + let xMax = -Infinity + let yMin = Infinity + let yMax = -Infinity + const record = ([x, y]: [number, number]) => { + if (!Number.isFinite(x) || !Number.isFinite(y)) return + xMin = Math.min(xMin, x) + xMax = Math.max(xMax, x) + yMin = Math.min(yMin, y) + yMax = Math.max(yMax, y) + } + for (let j = 0; j <= samples; j++) { + for (let i = 0; i <= samples; i++) { + record( + toSource( + box.x0 + (i / samples) * (box.x1 - box.x0), + box.y0 + (j / samples) * (box.y1 - box.y0) + ) + ) + } + } + for (const lat of [-90, 90]) { + const [px, py] = mapping.toWorld(0, lat) + if (px < box.x0 || px > box.x1 || py < box.y0 || py > box.y1) continue + // Every longitude meets at the pole. + for (const lng of [-180, -90, 0, 90, 180]) { + record(to4326.inverse(lng, lat)) + } + } + if (!Number.isFinite(xMin)) return null + const padX = (xMax - xMin) / samples + const padY = (yMax - yMin) / samples + return { + xMin: xMin - padX, + xMax: xMax + padX, + yMin: yMin - padY, + yMax: yMax + padY, + } +} + export function getRegionBounds({ regionX, regionY, @@ -394,17 +539,25 @@ export function selectLevelForZoom({ levels, projection, lonIndex, + worldFraction: mapWorldFraction, }: { mapZoom: number xyLimits: XYLimits | null levels: ResolutionLevel[] projection: ProjectionContext lonIndex?: number + /** + * The share of the map's world width the data spans, for a map whose + * world is not Web Mercator's. + */ + worldFraction?: number }): number { if (!xyLimits || levels.length === 0) return 0 const mapPixelsPerWorld = 256 * Math.pow(2, mapZoom) let worldFraction: number - if (projection.kind === 'epsg4326') { + if (mapWorldFraction !== undefined && mapWorldFraction > 0) { + worldFraction = mapWorldFraction + } else if (projection.kind === 'epsg4326') { // proj4's adjust_lon folds inputs outside [-180, 180] back into range // (e.g. forward(360,·)→0, forward(190,·)→-170), which collapses or // distorts the width for 0-360° stores and antimeridian-crossing diff --git a/src/region-renderer.render-path.test.ts b/src/region-renderer.render-path.test.ts index 7252775b..588c5633 100644 --- a/src/region-renderer.render-path.test.ts +++ b/src/region-renderer.render-path.test.ts @@ -13,6 +13,7 @@ import { maplibreFragmentShaderSource } from './shaders' import type { RenderContext } from './renderer-types' import type { ZarrRenderer } from './zarr-renderer' import type { RegionState } from './region-state' +import type { WorldMapping } from './world-mapping' /** * Per-frame render-path selection: the branch in RegionRenderer.render() that @@ -202,6 +203,12 @@ function mapboxContext( } } +/** Plate carrée over the top half of a registered projection's world square. */ +const PLATE: WorldMapping = { + toWorld: (lng, lat) => [(lng + 180) / 360, (90 - lat) / 360], + toLngLat: (x, y) => [x * 360 - 180, 90 - y * 360], +} + /** World offsets actually drawn, read back off the u_worldXOffset uploads. */ const drawnWorldOffsets = (gl: ReturnType) => gl.calls @@ -220,11 +227,17 @@ async function frame( { isGlobe = false, inactiveUniforms, - }: { isGlobe?: boolean; inactiveUniforms?: string[] } = {} + registeredProjection = false, + }: { + isGlobe?: boolean + inactiveUniforms?: string[] + registeredProjection?: boolean + } = {} ) { const renderer = await makeRenderer() seedRegion(renderer) renderer.onProjectionChange(isGlobe) + if (registeredProjection) renderer.setWorldMapping(PLATE) const gl = createRecordingGl({ inactiveUniforms }) const stub = stubRenderer(gl) @@ -299,6 +312,23 @@ describe('MapLibre render path', () => { expect(latLookup(gl)).toBe(1) }) + it('samples mesh texture coordinates in a registered projection', async () => { + const { gl } = await frame((g) => maplibreContext(g, 0), { + registeredProjection: true, + }) + expect(gl.callsTo('drawElements').length).toBeGreaterThan(0) + expect(latLookup(gl)).toBe(0) + }) + + it('drops the latitude lookup in a registered projection, tiles included, and restores it in Mercator', async () => { + const renderer = await makeRenderer() + const region = seedRegion(renderer) + renderer.setWorldMapping(PLATE) + expect(region.texLatEdges).toBe(null) + renderer.setWorldMapping(null) + expect(region.texLatEdges).toEqual([90, -90]) + }) + it('samples mesh texture coordinates on the globe', async () => { const { gl } = await frame((g) => maplibreContext(g, 0.5), { isGlobe: true, diff --git a/src/region-renderer.tile-lod.test.ts b/src/region-renderer.tile-lod.test.ts index f57dd6bb..67ed4f02 100644 --- a/src/region-renderer.tile-lod.test.ts +++ b/src/region-renderer.tile-lod.test.ts @@ -16,6 +16,8 @@ import { getRegionsInBounds, type RegionCoordinate } from './region-math' import type { ProjectionContext } from './projection-utils' import type { RegionState } from './region-state' import { boundsIntersect, getTileBounds } from './tile-utils' +import { latToMercatorNorm } from './map-utils' +import type { WorldMapping } from './world-mapping' /** * Per-tile level selection. The map picks a zoom per tile by distance when @@ -124,9 +126,22 @@ const failLevel = (renderer: RegionRenderer, levelIndex: number) => type Box = { x0: number; y0: number; x1: number; y1: number } +/** A box cut to a region's mesh: what of it the region can draw. */ +const withinMesh = (draw: Draw, box: Box): Box => { + const mesh = draw.region.meshBounds! + return { + x0: Math.max(box.x0, mesh.x0), + y0: Math.max(box.y0, mesh.y0), + x1: Math.min(box.x1, mesh.x1), + y1: Math.min(box.y1, mesh.y1), + } +} + /** * A draw's clip boxes back in absolute Mercator of the base world, for the - * passes drawn at one world offset. + * passes drawn at one world offset, each cut to the region's mesh: a box + * edge on the mesh edge is moved out, and only what lies inside the mesh + * is drawn. */ const absoluteBoxes = (draw: Draw, worldOffset = 0): Box[] => { const mesh = draw.region.meshBounds! @@ -137,12 +152,14 @@ const absoluteBoxes = (draw: Draw, worldOffset = 0): Box[] => { if (pass.worldOffset !== worldOffset) continue const p = pass.boxes for (let i = 0; i < p.length; i += 4) { - boxes.push({ - x0: p[i] + ax, - y0: p[i + 1] + ay, - x1: p[i + 2] + ax, - y1: p[i + 3] + ay, - }) + boxes.push( + withinMesh(draw, { + x0: p[i] + ax, + y0: p[i + 1] + ay, + x1: p[i + 2] + ax, + y1: p[i + 3] + ay, + }) + ) } } return boxes @@ -544,7 +561,9 @@ describe('per-tile level selection', () => { const tile = draw.region.levelIndex === 1 ? NEAR : FAR_EAST const boxes = absoluteBoxes(draw) expect(boxes).toHaveLength(1) - expect(sameBox(boxes[0], getTileBounds(tile.tileId))).toBe(true) + expect( + sameBox(boxes[0], withinMesh(draw, getTileBounds(tile.tileId))) + ).toBe(true) } }) @@ -568,12 +587,18 @@ describe('per-tile level selection', () => { if (draw.region.levelIndex === 0) { expect(offsets).toEqual([0]) expect( - sameBox(absoluteBoxes(draw, 0)[0], getTileBounds(FAR.tileId)) + sameBox( + absoluteBoxes(draw, 0)[0], + withinMesh(draw, getTileBounds(FAR.tileId)) + ) ).toBe(true) } else { expect(offsets).toEqual([-1]) expect( - sameBox(absoluteBoxes(draw, -1)[0], getTileBounds(NEAR.tileId)) + sameBox( + absoluteBoxes(draw, -1)[0], + withinMesh(draw, getTileBounds(NEAR.tileId)) + ) ).toBe(true) } } @@ -607,7 +632,9 @@ describe('per-tile level selection', () => { const tileBounds = getTileBounds(NEAR_IN_FAR_EAST.tileId) for (const draw of draws) { expect(cache.isProtected(draw.region.key)).toBe(true) - expect(sameBox(absoluteBoxes(draw)[0], tileBounds)).toBe(true) + expect( + sameBox(absoluteBoxes(draw)[0], withinMesh(draw, tileBounds)) + ).toBe(true) } await settle(renderer) }) @@ -1373,3 +1400,69 @@ describe('single-resolution datasets', () => { expect(wanted.has(makeRegionKey(0, 0, 0))).toBe(false) }) }) + +describe('registered projections', () => { + // Plate carrée over the top half of the world square, standing in for a + // map projection registered with addProjection. + const PLATE: WorldMapping = { + toWorld: (lng, lat) => [(lng + 180) / 360, (90 - lat) / 360], + toLngLat: (x, y) => { + const lat = 90 - y * 360 + return lat < -90 || lat > 90 ? null : [x * 360 - 180, lat] + }, + } + // The world square's tile 4/0/0: lon -180..-157.5, lat 90..67.5, two + // columns and four rows of the 32x32 grid. + const CORNER: TileRequest = { tileId: { z: 4, x: 0, y: 0 }, zoom: 0 } + + async function rendererInPlate() { + const made = await makeRenderer({ single: true }) + made.renderer.setWorldMapping(PLATE) + made.renderer.setCoveringTiles([CORNER]) + const map = mapAt(-180, 67.5, -157.5, 90, 0) + made.renderer.update(map, made.gl) + await settle(made.renderer) + return { ...made, map } + } + + it('selects the regions under each tile of the projection', async () => { + const { renderer } = await rendererInPlate() + + expect(wantedAt(renderer, 0)).toEqual([ + '0,0', + '0,1', + '0,2', + '0,3', + '1,0', + '1,1', + '1,2', + '1,3', + ]) + }) + + it('places regions in the world square of the projection', async () => { + const { renderer } = await rendererInPlate() + + // Region 1,3: lon -168.75..-157.5, lat 73.125..67.5. + const mesh = seam(renderer).regionCache.get(makeRegionKey(0, 1, 3))! + .meshBounds! + expect(mesh.x0).toBeCloseTo(11.25 / 360, 9) + expect(mesh.x1).toBeCloseTo(22.5 / 360, 9) + expect(mesh.y0).toBeCloseTo(16.875 / 360, 9) + expect(mesh.y1).toBeCloseTo(22.5 / 360, 9) + }) + + it('rebuilds cached regions in Web Mercator once the map leaves the projection', async () => { + const { renderer, invalidate } = await rendererInPlate() + const region = seam(renderer).regionCache.get(makeRegionKey(0, 1, 3))! + region.geometryUploaded = true + invalidate.mockClear() + + renderer.setWorldMapping(null) + + expect(region.geometryUploaded).toBe(false) + expect(region.meshBounds!.y1).toBeCloseTo(latToMercatorNorm(67.5), 6) + expect(region.mercatorBounds!.y1).toBeCloseTo(latToMercatorNorm(67.5), 6) + expect(invalidate).toHaveBeenCalled() + }) +}) diff --git a/src/region-renderer.ts b/src/region-renderer.ts index aef92249..41fb93da 100644 --- a/src/region-renderer.ts +++ b/src/region-renderer.ts @@ -62,11 +62,16 @@ import { getRegionBounds, getRegionSize, getRegionsInBounds, + getRegionsInWorldBox, selectLevelForZoom, type LonLatBounds, type RegionCoordinate, } from './region-math' -import { boundsIntersect, getTileLonLatBounds } from './tile-utils' +import { + boundsIntersect, + getTileBounds, + getTileLonLatBounds, +} from './tile-utils' import { tileClipBounds, tileClipPasses, type ClipTile } from './clip-boxes' import { RegionCache, @@ -88,7 +93,14 @@ import { import { RegionFetcher } from './region-fetcher' import { LevelLoader, type LevelLoadOutcome } from './level-loader' import { wrapError } from './errors' -import { createHybridMesh } from './mesh-reprojector' +import { createHybridMesh, createWorldMesh } from './mesh-reprojector' +import { + sourceToWorld, + worldBoundsOf, + worldFractionOf, + worldToSource, + type WorldMapping, +} from './world-mapping' import { type RequestCanceller, type LoadingManager, @@ -282,6 +294,15 @@ export class RegionRenderer { private rendersFromBandTextures: boolean = false // Track current projection for Mapbox's direct globe render path. private isGlobeProjection: boolean = false + // The map's world when it is in a projection registered with + // addProjection: regions are placed and selected in its world square, + // which its tiles divide, in place of normalized Web Mercator. + private worldMapping: WorldMapping | null = null + private worldFraction: { + mapping: WorldMapping + xyLimits: XYLimits + value: number | undefined + } | null = null // Fixed data scale for normalization (set at initialization, passed from ZarrLayer) private fixedDataScale: number = 1 // The layer's `fillValue` option, which replaces the store's fill. @@ -682,6 +703,18 @@ export class RegionRenderer { const memo = this.tileRegionsMemo.get(key) if (memo) return memo if (this.tileRegionsMemo.size > 4096) this.tileRegionsMemo.clear() + if (this.worldMapping) { + const regions = getRegionsInWorldBox({ + box: getTileBounds(tileId), + xyLimits: this.xyLimits, + levelMeta, + projection: this.projection, + mapping: this.worldMapping, + latIsAscending: this.latIsAscending, + }) + this.tileRegionsMemo.set(key, regions) + return regions + } const regions: RegionCoordinate[] = [] const seen = new Set() for (const bounds of extents ?? [ @@ -932,9 +965,7 @@ export class RegionRenderer { }) } - const { lng, lat } = this.mapCenter(map) - const cx = (lng + 180) / 360 - const cy = latToMercatorNorm(lat) + const [cx, cy] = this.centerInWorld(map) const entries: Array<{ tile: TileRequest key: string @@ -1030,7 +1061,20 @@ export class RegionRenderer { } } - /** Regions ordered by their distance from the map center, in degrees. */ + /** + * The map center in the world its tiles divide: the registered + * projection's world square, else normalized Web Mercator. + */ + private centerInWorld(map: MapLike): [number, number] { + const { lng, lat } = this.mapCenter(map) + if (this.worldMapping) return this.worldMapping.toWorld(lng, lat) + return [(lng + 180) / 360, latToMercatorNorm(lat)] + } + + /** + * Regions ordered by their distance from the map center: in degrees, or + * in the world square of a registered projection. + */ private nearestCenterFirst( regions: RegionCoordinate[], runtime: LevelRuntime, @@ -1038,10 +1082,18 @@ export class RegionRenderer { ): RegionCoordinate[] { const { lng, lat } = this.mapCenter(map) const inverse = this.projection.toWGS84?.inverse + const mapping = this.worldMapping + const to4326 = this.projection.to4326 + const [cx, cy] = this.centerInWorld(map) const distance = (region: RegionCoordinate) => { const b = this.getRegionBounds(region.regionX, region.regionY, runtime) const xMid = (b.xMin + b.xMax) / 2 const yMid = (b.yMin + b.yMax) / 2 + if (mapping && to4326) { + const [x, y] = sourceToWorld(mapping, to4326)(xMid, yMid) + const d = (x - cx) ** 2 + (y - cy) ** 2 + return Number.isFinite(d) ? d : Infinity + } const [lon, latitude] = inverse ? inverse(xMid, yMid) : [xMid, yMid] let dLon = ((lon - lng + 540) % 360) - 180 if (!Number.isFinite(dLon)) dLon = 180 @@ -1496,30 +1548,44 @@ export class RegionRenderer { const lonSubdivisions = subdivisionsForSpan(lonSpan) const latSubdivisions = subdivisionsForSpan(latSpan) - const meshResult = createHybridMesh({ - geoBounds, - width: region.width, - height: region.height, - lonSubdivisions, - latSubdivisions, - transformer: this.projection.to4326, - latIsAscending: this.latIsAscending, - allowUnwrappedLongitudes: this.projection.kind === 'epsg4326', - // EPSG:4326 latitude is resolved per fragment on the flat map - // (texLatEdges), so only other projections refine against Mercator. - refineInMercator: this.projection.kind !== 'epsg4326', - // The mesh encodes vertices as deltas from a per-region mercator origin - // (deriveLocalMercAnchor); renderRegion uploads a matching per-region - // anchor_clip, keeping the eye origin near the on-screen region for - // high-zoom precision (no pan/zoom jitter, no seams). See - // VERTEX_TO_WGS84_TO_MERCATOR. - }) + const meshResult = this.worldMapping + ? createWorldMesh({ + geoBounds, + width: region.width, + height: region.height, + lonSubdivisions, + latSubdivisions, + toWorld: sourceToWorld(this.worldMapping, this.projection.to4326), + fromWorld: worldToSource(this.worldMapping, this.projection.to4326), + latIsAscending: this.latIsAscending, + allowUnwrappedLongitudes: this.projection.kind === 'epsg4326', + }) + : createHybridMesh({ + geoBounds, + width: region.width, + height: region.height, + lonSubdivisions, + latSubdivisions, + transformer: this.projection.to4326, + latIsAscending: this.latIsAscending, + allowUnwrappedLongitudes: this.projection.kind === 'epsg4326', + // EPSG:4326 latitude is resolved per fragment on the flat map + // (texLatEdges), so only other projections refine against Mercator. + refineInMercator: this.projection.kind !== 'epsg4326', + // The mesh encodes vertices as deltas from a per-region mercator origin + // (deriveLocalMercAnchor); renderRegion uploads a matching per-region + // anchor_clip, keeping the eye origin near the on-screen region for + // high-zoom precision (no pan/zoom jitter, no seams). See + // VERTEX_TO_WGS84_TO_MERCATOR. + }) region.vertexArr = meshResult.positions region.pixCoordArr = meshResult.texCoords region.indexArr = meshResult.indices region.meshBounds = meshResult.meshBounds region.indexCount = region.indexArr.length - if (this.projection.kind === 'epsg4326') { + // The latitude lookup reads positions as Mercator, which a registered + // projection's world square is not. + if (this.projection.kind === 'epsg4326' && !this.worldMapping) { region.texLatEdges = this.latIsAscending ? [geoBounds.yMin, geoBounds.yMax] : [geoBounds.yMax, geoBounds.yMin] @@ -1818,9 +1884,38 @@ export class RegionRenderer { levels: this.levels, projection: this.projection, lonIndex: this.dimIndices.lon?.index, + worldFraction: this.dataWorldFraction(), }) } + /** + * The share of a registered projection's world width the data spans, or + * undefined on a Web Mercator world. + */ + private dataWorldFraction(): number | undefined { + if (!this.worldMapping || !this.xyLimits || !this.projection.to4326) { + return undefined + } + const cached = this.worldFraction + if ( + cached?.mapping === this.worldMapping && + cached.xyLimits === this.xyLimits + ) { + return cached.value + } + const value = worldFractionOf( + this.xyLimits, + this.worldMapping, + this.projection.to4326 + ) + this.worldFraction = { + mapping: this.worldMapping, + xyLimits: this.xyLimits, + value, + } + return value + } + render(renderer: ZarrRenderer, context: RenderContext): void { context = { ...context, @@ -2069,15 +2164,49 @@ export class RegionRenderer { return computeMercatorBoundsFromProjection(this.xyLimits, this.projection) } + /** + * A region's box in the world the map's tiles divide: the registered + * projection's world square, else normalized Web Mercator. + */ private computeRegionMercatorBounds(bounds: { xMin: number xMax: number yMin: number yMax: number }): MercatorBounds { + if (this.worldMapping && this.projection.to4326) { + return ( + worldBoundsOf(bounds, this.worldMapping, this.projection.to4326) ?? { + x0: 0, + y0: 0, + x1: 0, + y1: 0, + } + ) + } return computeRegionMercatorBounds(bounds, this.projection) } + /** + * Place regions in the world of a map in a registered projection, or in + * Web Mercator for null. Cached regions are rebuilt for the new world, + * keeping their data, and the tile plans start over. + */ + setWorldMapping(mapping: WorldMapping | null): void { + if (mapping === this.worldMapping) return + this.worldMapping = mapping + this.tileRegionsMemo.clear() + this.tileLevels = new Map() + this.lastViewportHash = '' + for (const region of this.regionCache.values()) { + if (!region.levelMeta) continue + if (!region.mercatorBounds && !region.vertexArr) continue + region.mercatorBounds = null + this.createRegionGeometry(region.regionX, region.regionY, region) + } + this.invalidate() + } + async setSelector(selector: NormalizedSelector): Promise { this.selector = selector this.bandNames = getBands(this.variable, selector) diff --git a/src/types.ts b/src/types.ts index f5a2a6c5..6190b07f 100644 --- a/src/types.ts +++ b/src/types.ts @@ -1,6 +1,7 @@ import * as zarr from 'zarrita' import type { XYLimits } from './map-utils' import type { UniformValue } from './renderer-types' +import type { MapProjection } from './world-mapping' /** Bounds tuple: [xMin, yMin, xMax, yMax] */ export type Bounds = [number, number, number, number] @@ -175,6 +176,15 @@ export interface ZarrLayerOptions { * Default: `false` */ renderPoles?: boolean + /** + * Projections registered with MapLibre's experimental `addProjection` + * (maplibre-gl >= 6.15) that the layer can render in: pass the same + * objects given to `addProjection`. While the map is in a registered + * projection missing from this list, the layer is hidden and warns once. + * + * Default: `[]` + */ + mapProjections?: MapProjection[] } export type CRS = 'EPSG:4326' | 'EPSG:3857' diff --git a/src/world-mapping.test.ts b/src/world-mapping.test.ts new file mode 100644 index 00000000..cf95a367 --- /dev/null +++ b/src/world-mapping.test.ts @@ -0,0 +1,741 @@ +import { describe, it, expect } from 'vitest' +import proj4 from 'proj4' +import { + sourceToWorld, + worldBoundsOf, + worldMappingFor, + worldToSource, + worldFractionOf, + type MapProjection, + type WorldMapping, +} from './world-mapping' +import { getRegionsInWorldBox, selectLevelForZoom } from './region-math' +import { createWorldMesh } from './mesh-reprojector' +import { createProjectionContext } from './projection-utils' +import type { MapLike } from './types' + +/** + * A map in a projection registered with MapLibre's `addProjection` places + * custom layers in the world square of its tile matrix set. These fixtures + * stand in for a registered projection's mapping with closed forms. + */ + +/** Plate carrée over the top half of the world square. */ +const PLATE: WorldMapping = { + toWorld: (lng, lat) => [(lng + 180) / 360, (90 - lat) / 360], + toLngLat: (x, y) => { + const lat = 90 - y * 360 + return lat < -90 || lat > 90 ? null : [x * 360 - 180, lat] + }, +} + +/** A south polar azimuthal projection, centered in the world square. */ +const POLAR: WorldMapping = { + toWorld: (lng, lat) => { + const r = (lat + 90) / 360 + const a = (lng * Math.PI) / 180 + return [0.5 + r * Math.sin(a), 0.5 - r * Math.cos(a)] + }, + toLngLat: (x, y) => { + const dx = x - 0.5 + const dy = 0.5 - y + return [(Math.atan2(dx, dy) * 180) / Math.PI, Math.hypot(dx, dy) * 360 - 90] + }, +} + +/** + * South polar stereographic: the north pole lands a finite but enormous + * distance away, as proj4 places it. + */ +const STEREO: WorldMapping = { + toWorld: (lng, lat) => { + const r = Math.tan((((lat + 90) / 2) * Math.PI) / 180) / 4 + const a = (lng * Math.PI) / 180 + return [0.5 + r * Math.sin(a), 0.5 - r * Math.cos(a)] + }, + toLngLat: (x, y) => { + const dx = x - 0.5 + const dy = 0.5 - y + return [ + (Math.atan2(dx, dy) * 180) / Math.PI, + (Math.atan(Math.hypot(dx, dy) * 4) * 360) / Math.PI - 90, + ] + }, +} + +const IDENTITY_4326 = { + forward: (x: number, y: number): [number, number] => [x, y], + inverse: (x: number, y: number): [number, number] => [x, y], +} + +/** Plate carrée in degrees, tile 0/0/0 spanning lon -180..180 from lat 90 down. */ +const plateCarree = ( + overrides: Partial = {} +): MapProjection => ({ + name: 'custom:plate', + projection: { + forward: ([lng, lat]) => [lng, lat], + inverse: ([x, y]) => [x, y], + ...overrides, + }, + tileMatrixSet: { origin: [-180, 90], extentAtZoom0: 360 }, +}) + +const mapIn = (type: string) => + ({ getProjection: () => ({ type }) } as unknown as MapLike) + +describe('worldMappingFor', () => { + it('maps the projection the map is in from its definition', () => { + const mapping = worldMappingFor(mapIn('custom:plate'), [plateCarree()])! + + // World y grows downwards from the tile matrix set's top-left origin. + expect(mapping.toWorld(0, 90)).toEqual([0.5, 0]) + expect(mapping.toWorld(-180, 0)).toEqual([0, 0.25]) + expect(mapping.toLngLat(0.75, 0.125)).toEqual([90, 45]) + }) + + it('gives the same mapping each frame for the same projection', () => { + const projections = [plateCarree()] + expect(worldMappingFor(mapIn('custom:plate'), projections)).toBe( + worldMappingFor(mapIn('custom:plate'), projections) + ) + }) + + it('returns null for built-in projections and ones it was not given', () => { + const projections = [plateCarree()] + expect(worldMappingFor(mapIn('mercator'), projections)).toBe(null) + expect(worldMappingFor(mapIn('globe'), projections)).toBe(null) + expect(worldMappingFor(mapIn('EPSG:3031'), projections)).toBe(null) + expect(worldMappingFor(null, projections)).toBe(null) + }) + + it('reports positions the projection cannot place', () => { + const mapping = worldMappingFor(mapIn('custom:plate'), [ + plateCarree({ + forward: ([lng, lat]) => { + if (lng > 90) throw new Error('outside') + return [lng > 45 ? Infinity : lng, lat] + }, + inverse: ([x, y]) => (x > 0 ? null : [x, y]), + }), + ])! + + expect(mapping.toWorld(60, 0)).toEqual([NaN, NaN]) + expect(mapping.toWorld(120, 0)).toEqual([NaN, NaN]) + expect(mapping.toLngLat(0.75, 0.25)).toBe(null) + }) + + it('has no answer where the inverse does not lead back to the position', () => { + const mapping = worldMappingFor(mapIn('custom:plate'), [ + plateCarree({ inverse: ([x, y]) => (x > 0 ? [x - 360, y] : [x, y]) }), + ])! + + expect(mapping.toLngLat(0.25, 0.25)).toEqual([-90, 0]) + expect(mapping.toLngLat(0.75, 0.25)).toBe(null) + }) + + it('takes an inverse that only approximately leads back, as MapLibre does', () => { + const rounded = (v: number) => Math.round(v * 100) / 100 + const mapping = worldMappingFor(mapIn('custom:plate'), [ + plateCarree({ inverse: ([x, y]) => [rounded(x), rounded(y)] }), + ])! + + expect(mapping.toLngLat(0.5 + 0.004 / 360, 0.25)).toEqual([0, 0]) + }) +}) + +describe('worldBoundsOf', () => { + it('maps a rectangle under a linear mapping exactly', () => { + expect( + worldBoundsOf( + { xMin: -180, xMax: 0, yMin: 0, yMax: 90 }, + PLATE, + IDENTITY_4326 + ) + ).toEqual({ x0: 0, y0: 0, x1: 0.5, y1: 0.25 }) + }) + + it('reaches the curved edges of a cap around the pole', () => { + const bounds = worldBoundsOf( + { xMin: -180, xMax: 180, yMin: -90, yMax: -72 }, + POLAR, + IDENTITY_4326 + )! + const r = 18 / 360 + + for (const [edge, cap] of [ + [0.5 - bounds.x0, r], + [bounds.x1 - 0.5, r], + [0.5 - bounds.y0, r], + [bounds.y1 - 0.5, r], + ]) { + expect(edge).toBeGreaterThanOrEqual(cap) + expect(edge).toBeLessThan(cap * 1.05) + } + }) + + it('holds a curved edge that bulges out between samples', () => { + // The equator from lon 10 to 190 reaches farthest east at lon 90, between + // the samples at 55 and 100. + const bounds = worldBoundsOf( + { xMin: 10, xMax: 190, yMin: -90, yMax: 0 }, + POLAR, + IDENTITY_4326, + 4 + )! + + expect(bounds.x1).toBeGreaterThanOrEqual(0.75) + expect(bounds.x1).toBeLessThan(0.76) + }) + + it('is null when nothing lands in the world', () => { + const nowhere: WorldMapping = { ...PLATE, toWorld: () => [NaN, NaN] } + expect( + worldBoundsOf( + { xMin: 0, xMax: 1, yMin: 0, yMax: 1 }, + nowhere, + IDENTITY_4326 + ) + ).toBe(null) + }) +}) + +describe('getRegionsInWorldBox', () => { + const projection = createProjectionContext({ + crs: 'EPSG:4326', + proj4def: null, + xyLimits: { xMin: -180, xMax: 180, yMin: -90, yMax: 90 }, + }) + // 4x4 regions of 90 by 45 degrees. + const levelMeta = { + width: 32, + height: 16, + regionSize: [4, 8] as [number, number], + } + const regionsIn = ( + box: { x0: number; y0: number; x1: number; y1: number }, + mapping = PLATE + ) => + getRegionsInWorldBox({ + box, + xyLimits: { xMin: -180, xMax: 180, yMin: -90, yMax: 90 }, + levelMeta, + projection, + mapping, + latIsAscending: false, + }) + .map(({ regionX, regionY }) => `${regionX},${regionY}`) + .sort() + + it('selects the regions under a tile, leaving out those that only touch it', () => { + // Tile 1/0/0 is the western hemisphere in this world. + expect(regionsIn({ x0: 0, y0: 0, x1: 0.5, y1: 0.5 })).toEqual([ + '0,0', + '0,1', + '0,2', + '0,3', + '1,0', + '1,1', + '1,2', + '1,3', + ]) + }) + + it('keeps a region whose curved edge bulges into the box between samples', () => { + // One region, lon 10..190 and lat -90..0; its equator edge reaches past + // x 0.748 only around lon 90. + const limits = { xMin: 10, xMax: 190, yMin: -90, yMax: 0 } + expect( + getRegionsInWorldBox({ + box: { x0: 0.747, y0: 0.499, x1: 0.749, y1: 0.501 }, + xyLimits: limits, + levelMeta: { width: 1, height: 1, regionSize: [1, 1] }, + projection: createProjectionContext({ + crs: 'EPSG:4326', + proj4def: null, + xyLimits: limits, + }), + mapping: POLAR, + latIsAscending: false, + }) + ).toEqual([{ regionX: 0, regionY: 0 }]) + }) + + it('selects regions of a 0-360° dataset west of the antimeridian', () => { + const limits = { xMin: 0, xMax: 360, yMin: -90, yMax: 90 } + // Lon 260 (-100) at 60°S, in region 2,3 (lon 180..270, lat -45..-90). + const [x, y] = POLAR.toWorld(-100, -60) + expect( + getRegionsInWorldBox({ + box: { x0: x - 0.002, y0: y - 0.002, x1: x + 0.002, y1: y + 0.002 }, + xyLimits: limits, + levelMeta, + projection: createProjectionContext({ + crs: 'EPSG:4326', + proj4def: null, + xyLimits: limits, + }), + mapping: POLAR, + latIsAscending: false, + }) + ).toEqual([{ regionX: 2, regionY: 3 }]) + }) + + it('selects the ring of regions around the pole of a polar world', () => { + // A box at the pole holds the southernmost row, every longitude. + expect( + regionsIn({ x0: 0.45, y0: 0.45, x1: 0.55, y1: 0.55 }, POLAR) + ).toEqual(['0,3', '1,3', '2,3', '3,3']) + }) +}) + +describe('worldFractionOf', () => { + it('is the share of the world width the data spans', () => { + expect( + worldFractionOf( + { xMin: -180, xMax: 0, yMin: 0, yMax: 90 }, + PLATE, + IDENTITY_4326 + ) + ).toBeCloseTo(0.5, 9) + }) + + it('counts only the world square, past which far points are placed', () => { + // Global data reaches the north pole, which the projection places far + // outside the world square. + expect( + worldFractionOf( + { xMin: -180, xMax: 180, yMin: -90, yMax: 90 }, + STEREO, + IDENTITY_4326 + ) + ).toBeCloseTo(1, 1) + }) + + it('counts only the part of a narrow strip inside the world square', () => { + // A one-degree strip reaching the north pole: inside the square it is a + // sliver from the south pole out to the square's edge. + const fraction = worldFractionOf( + { xMin: 0, xMax: 1, yMin: -90, yMax: 90 }, + STEREO, + IDENTITY_4326 + )! + + expect(fraction).toBeGreaterThan(0.005) + expect(fraction).toBeLessThanOrEqual(0.5 * Math.sin(Math.PI / 180)) + }) + + it('is undefined when nothing lands in the world', () => { + expect( + worldFractionOf( + { xMin: 0, xMax: 1, yMin: 0, yMax: 1 }, + { ...PLATE, toWorld: () => [NaN, NaN] }, + IDENTITY_4326 + ) + ).toBe(undefined) + }) +}) + +describe('selectLevelForZoom in a registered projection', () => { + const projection = createProjectionContext({ + crs: 'EPSG:4326', + proj4def: null, + xyLimits: { xMin: -180, xMax: 180, yMin: -90, yMax: 90 }, + }) + const levels = [ + { asset: '0', shape: [128, 256] }, + { asset: '1', shape: [256, 512] }, + ] + const select = (worldFraction?: number) => + selectLevelForZoom({ + mapZoom: 1, + xyLimits: { xMin: -180, xMax: 180, yMin: -90, yMax: 90 }, + levels, + projection, + worldFraction, + }) + + it('reads coarser data when the data spans less of the world', () => { + expect(select()).toBe(1) + expect(select(1)).toBe(1) + expect(select(0.5)).toBe(0) + }) +}) + +describe('createWorldMesh', () => { + const geoBounds = { xMin: -180, xMax: 0, yMin: 0, yMax: 90 } + const mesh = (mapping: WorldMapping) => + createWorldMesh({ + geoBounds, + width: 16, + height: 8, + lonSubdivisions: 4, + latSubdivisions: 4, + toWorld: (x, y) => mapping.toWorld(x, y), + fromWorld: (x, y) => mapping.toLngLat(x, y) ?? [NaN, NaN], + latIsAscending: false, + }) + const worldPositions = (result: ReturnType) => { + const { x0, y0, x1, y1 } = result.meshBounds + const points: Array<[number, number]> = [] + for (let i = 0; i < result.positions.length; i += 2) { + points.push([ + (x0 + x1) / 2 + (result.positions[i] * (x1 - x0)) / 2, + (y0 + y1) / 2 + (result.positions[i + 1] * (y1 - y0)) / 2, + ]) + } + return points + } + + it('places vertices at their world positions, around the region anchor', () => { + const result = mesh(PLATE) + + expect(result.meshBounds.x0).toBeCloseTo(0, 9) + expect(result.meshBounds.x1).toBeCloseTo(0.5, 9) + expect(result.meshBounds.y0).toBeCloseTo(0, 9) + expect(result.meshBounds.y1).toBeCloseTo(0.25, 9) + const points = worldPositions(result) + for (let i = 0; i < points.length; i++) { + const u = result.texCoords[i * 2] + const v = result.texCoords[i * 2 + 1] + // v runs north to south, as the data rows do. + expect(points[i][0]).toBeCloseTo(u * 0.5, 6) + expect(points[i][1]).toBeCloseTo(v * 0.25, 6) + } + expect(result.indices.length).toBeGreaterThan(0) + }) + + it('drops triangles with a vertex the projection cannot place', () => { + const western: WorldMapping = { + ...PLATE, + toWorld: (lng, lat) => (lng > -90 ? [NaN, NaN] : PLATE.toWorld(lng, lat)), + } + const result = mesh(western) + const points = worldPositions(result) + + expect(result.indices.length).toBeGreaterThan(0) + for (const index of result.indices) { + expect(points[index][0]).toBeLessThanOrEqual(0.25 + 1e-9) + } + expect(result.meshBounds.x1).toBeCloseTo(0.25, 9) + }) + + it('meets the neighboring region along a curved edge without a gap, at any density', () => { + const region = (yMin: number, yMax: number, lonSubdivisions: number) => + createWorldMesh({ + geoBounds: { xMin: -90, xMax: 0, yMin, yMax }, + width: 16, + height: 8, + lonSubdivisions, + latSubdivisions: 4, + toWorld: (x, y) => POLAR.toWorld(x, y), + fromWorld: (x, y) => POLAR.toLngLat(x, y) ?? [NaN, NaN], + latIsAscending: true, + }) + // The edge at 60°S, a circle around the pole. + const edge = (result: ReturnType, v: number) => { + const points = worldPositions(result) + return [...new Set(result.indices)] + .filter((i) => result.texCoords[i * 2 + 1] === v) + .map((i) => points[i]) + } + const onPolyline = ( + [x, y]: [number, number], + polyline: Array<[number, number]> + ) => + polyline.some(([ax, ay], k) => { + const [bx, by] = polyline[(k + 1) % polyline.length] + const t = Math.max( + 0, + Math.min( + 1, + ((x - ax) * (bx - ax) + (y - ay) * (by - ay)) / + ((bx - ax) ** 2 + (by - ay) ** 2 || 1) + ) + ) + return Math.hypot(x - ax - t * (bx - ax), y - ay - t * (by - ay)) < 1e-6 + }) + const byAngle = (points: Array<[number, number]>) => + [...points].sort( + (a, b) => + Math.atan2(a[0] - 0.5, 0.5 - a[1]) - + Math.atan2(b[0] - 0.5, 0.5 - b[1]) + ) + + // Neighbors meshed at the same density and at different ones. + for (const [innerCells, outerCells] of [ + [4, 4], + [4, 7], + ]) { + const inner = edge(region(-90, -60, innerCells), 1) + const outer = edge(region(-60, -30, outerCells), 0) + for (const point of outer) + expect(onPolyline(point, byAngle(inner))).toBe(true) + for (const point of inner) + expect(onPolyline(point, byAngle(outer))).toBe(true) + } + }) + + it('meshes a region reaching the pole the projection sends far away, past the grid', () => { + // Antarctic polar stereographic places the north pole ~1e23 m out, where + // proj4's inverse no longer leads back. + const antarctic = worldMappingFor(mapIn('EPSG:3031'), [ + { + 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, + }, + }, + ])! + const result = createWorldMesh({ + geoBounds: { xMin: -36, xMax: 36, yMin: 30, yMax: 90 }, + width: 360, + height: 300, + lonSubdivisions: 8, + latSubdivisions: 8, + toWorld: (x, y) => antarctic.toWorld(x, y), + fromWorld: (x, y) => antarctic.toLngLat(x, y) ?? [NaN, NaN], + latIsAscending: true, + }) + + // All of it lies past the edge of the tile matrix set. + expect(result.indices.length).toBe(0) + }) + + it('keeps the part of a region inside a grid far smaller than the region', () => { + // A regional grid: the world square spans 36 by 36 degrees, so the + // western hemisphere covers five squares across and its cells are wider + // than the square. + const regional: WorldMapping = { + toWorld: (lng, lat) => [(lng + 108) / 36, (54 - lat) / 36], + toLngLat: (x, y) => [x * 36 - 108, 54 - y * 36], + } + const result = mesh(regional) + const points = worldPositions(result) + + expect(result.indices.length).toBeGreaterThan(0) + for (let t = 0; t < result.indices.length; t += 3) { + const xs = [0, 1, 2].map((k) => points[result.indices[t + k]][0]) + const ys = [0, 1, 2].map((k) => points[result.indices[t + k]][1]) + expect(Math.max(...xs)).toBeGreaterThanOrEqual(0) + expect(Math.min(...xs)).toBeLessThanOrEqual(1) + expect(Math.max(...ys)).toBeGreaterThanOrEqual(0) + expect(Math.min(...ys)).toBeLessThanOrEqual(1) + } + }) + + it('drops triangles crossing a seam, however short their edges', () => { + // Lambert conformal conic cuts the world along lon 180. Near the cone's + // tip, triangles crossing the cut have edges shorter than most others. + const lcc = proj4( + 'EPSG:4326', + '+proj=lcc +lat_0=0 +lon_0=0 +lat_1=30 +lat_2=60 +datum=WGS84 +units=m' + ) + const toWorld = (lng: number, lat: number): [number, number] => { + const [x, y] = lcc.forward([lng, lat]) + return [0.5 + x / 4e7, 0.75 - y / 4e7] + } + const result = createWorldMesh({ + geoBounds: { xMin: 0, xMax: 360, yMin: 0, yMax: 90 }, + width: 360, + height: 90, + lonSubdivisions: 128, + latSubdivisions: 90, + toWorld, + fromWorld: (x, y) => { + const [lng, lat] = lcc.inverse([(x - 0.5) * 4e7, (0.75 - y) * 4e7]) + const [bx, by] = toWorld(lng, lat) + return Math.hypot(bx - x, by - y) <= 1e-3 ? [lng, lat] : [NaN, NaN] + }, + latIsAscending: true, + allowUnwrappedLongitudes: true, + }) + + expect(result.indices.length).toBeGreaterThan(0) + for (let t = 0; t < result.indices.length; t += 3) { + const lngs = [0, 1, 2].map( + (k) => result.texCoords[result.indices[t + k] * 2] * 360 + ) + expect(Math.min(...lngs) <= 180 && Math.max(...lngs) > 180).toBe(false) + } + }) + + // Equal Earth bends hard near its pole line; a cone with standard parallels + // near the equator stretches without bound toward its tip. Long diagonal + // edges there stray from their chords without crossing anything. + it.each([ + ['Equal Earth', '+proj=eqearth +datum=WGS84', 90, 180, 34487920], + [ + 'a sharp cone', + '+proj=lcc +lat_0=90 +lon_0=0 +lat_1=5 +lat_2=10 +datum=WGS84', + 0, + 90, + 1e8, + ], + ])( + 'keeps triangles where %s bends sharply but has no seam', + (_, definition, lngMin, lngMax, scale) => { + const converter = proj4('EPSG:4326', definition) + const result = createWorldMesh({ + geoBounds: { xMin: lngMin, xMax: lngMax, yMin: 60, yMax: 90 }, + width: 256, + height: 256, + lonSubdivisions: 4, + latSubdivisions: 4, + toWorld: (lng, lat) => { + const [x, y] = converter.forward([lng, lat]) + return [0.5 + x / scale, 0.5 - y / scale] + }, + fromWorld: (x, y) => + converter.inverse([(x - 0.5) * scale, (0.5 - y) * scale]), + latIsAscending: true, + }) + + let area = 0 + for (let t = 0; t < result.indices.length; t += 3) { + const [a, b, c] = [0, 1, 2].map((k) => result.indices[t + k] * 2) + const uv = result.texCoords + area += + Math.abs( + (uv[b] - uv[a]) * (uv[c + 1] - uv[a + 1]) - + (uv[c] - uv[a]) * (uv[b + 1] - uv[a + 1]) + ) / 2 + } + expect(area).toBeCloseTo(1, 9) + } + ) + + describe('data the map cannot show', () => { + const SOUTH_POLAR = + '+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' + const ARCTIC = + '+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' + 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' + const registered = ( + definition: string, + origin: [number, number], + extentAtZoom0: number + ) => + worldMappingFor(mapIn('custom:test'), [ + { + name: 'custom:test', + projection: proj4('EPSG:4326', definition), + tileMatrixSet: { origin, extentAtZoom0 }, + }, + ])! + const meshOf = ( + mapping: WorldMapping, + sourceDef: string, + [xMin, yMin, xMax, yMax]: number[] + ) => { + const { to4326 } = createProjectionContext({ + crs: 'EPSG:4326', + proj4def: sourceDef, + xyLimits: { xMin, xMax, yMin, yMax }, + }) + return createWorldMesh({ + geoBounds: { xMin, xMax, yMin, yMax }, + width: 128, + height: 128, + lonSubdivisions: 8, + latSubdivisions: 8, + toWorld: sourceToWorld(mapping, to4326!), + fromWorld: worldToSource(mapping, to4326!), + latIsAscending: false, + }) + } + const overlapsWorld = (result: ReturnType) => { + const points = worldPositions(result) + for (let t = 0; t < result.indices.length; t += 3) { + const p = [0, 1, 2].map((k) => points[result.indices[t + k]]) + const xs = p.map(([x]) => x) + const ys = p.map(([, y]) => y) + if ( + Math.max(...xs) >= 0 && + Math.min(...xs) <= 1 && + Math.max(...ys) >= 0 && + Math.min(...ys) <= 1 + ) { + return true + } + } + return false + } + // Antarctic data around the south pole, in polar stereographic metres. + const ANTARCTICA = [-3315363, -3316901, 3316909, 3315371] + + it.each([ + // The south pole lands at infinity; triangles around it would wrap + // the whole map. + ['an Arctic map', ARCTIC, [-4194304, 4194304], 8388608], + ['a CONUS map', CONUS_ALBERS, [-3000000, 4750000], 6000000], + ] as const)( + 'draws nothing of Antarctica on %s', + (_, definition, origin, extent) => { + const mapping = registered(definition, [...origin], extent) + expect(overlapsWorld(meshOf(mapping, SOUTH_POLAR, ANTARCTICA))).toBe( + false + ) + } + ) + + it('keeps a triangle crossing the map with every corner off it', () => { + // Coarse global lon/lat data in CONUS Albers: no vertex falls inside + // the map, so the triangles spanning it are what draws it. + const mapping = registered(CONUS_ALBERS, [-3000000, 4750000], 6000000) + const result = createWorldMesh({ + geoBounds: { xMin: -180, xMax: 180, yMin: -90, yMax: 90 }, + width: 128, + height: 64, + lonSubdivisions: 2, + latSubdivisions: 2, + toWorld: (lng, lat) => mapping.toWorld(lng, lat), + fromWorld: (x, y) => mapping.toLngLat(x, y) ?? [NaN, NaN], + latIsAscending: false, + allowUnwrappedLongitudes: true, + maxError: Infinity, + }) + const points = worldPositions(result) + const offMap = ([x, y]: [number, number]) => + x < 0 || x > 1 || y < 0 || y > 1 + const coversCenter = ([a, b, c]: [number, number][]) => { + const side = (p: [number, number], q: [number, number]) => + (q[0] - p[0]) * (0.5 - p[1]) - (q[1] - p[1]) * (0.5 - p[0]) + const s = [side(a, b), side(b, c), side(c, a)] + return s.every((v) => v >= 0) || s.every((v) => v <= 0) + } + const covering = [] + for (let t = 0; t < result.indices.length; t += 3) { + const corners = [0, 1, 2].map((k) => points[result.indices[t + k]]) + if (coversCenter(corners)) covering.push(corners) + } + + expect(covering.length).toBeGreaterThan(0) + expect(covering.every((corners) => corners.every(offMap))).toBe(true) + }) + }) + + it('drops triangles spanning a seam in the projection', () => { + // East of lon -90 the world wraps to the far side, as past a seam. + const seam: WorldMapping = { + ...PLATE, + toWorld: (lng, lat) => { + const [x, y] = PLATE.toWorld(lng, lat) + return [lng > -90 ? x + 1 : x, y] + }, + } + const result = mesh(seam) + const points = worldPositions(result) + + for (let t = 0; t < result.indices.length; t += 3) { + const xs = [0, 1, 2].map((k) => points[result.indices[t + k]][0]) + expect(Math.max(...xs) - Math.min(...xs)).toBeLessThan(0.5) + } + }) +}) diff --git a/src/world-mapping.ts b/src/world-mapping.ts new file mode 100644 index 00000000..841401c8 --- /dev/null +++ b/src/world-mapping.ts @@ -0,0 +1,242 @@ +/** + * @module world-mapping + * + * The flat world of a map in a projection registered with MapLibre's + * `addProjection` (maplibre-gl >= 6.15). Such a map lays its tiles over a + * square tile matrix set in the CRS's own plane, and custom layers take + * positions in that square: [0, 0] is the top-left of tile 0/0/0 and [1, 1] + * its bottom-right, with y growing downwards like tile rows. The layer + * places regions there in place of normalized Web Mercator, and the tile + * math it shares with Mercator (tile ids, clip boxes) carries over unchanged. + */ + +import { customProjectionName } from './map-utils' +import type { MapLike } from './types' + +/** + * A projection registered with MapLibre's `addProjection`, in the shape of + * its `CRS` type, so the same object can be passed to both. + */ +export interface MapProjection { + /** The registered name, e.g. `EPSG:3031`. */ + name: string + /** WGS84 degrees to CRS units and back, like a proj4 converter. */ + projection: { + forward(position: number[]): number[] + inverse(position: number[]): number[] | null + } + tileMatrixSet: { + /** Top-left corner of tile 0/0/0, in CRS units. */ + origin: [number, number] + /** Width of tile 0/0/0, in CRS units. */ + extentAtZoom0: number + } +} + +export interface WorldMapping { + /** WGS84 degrees to the world square; NaN where the CRS has no answer. */ + toWorld(lng: number, lat: number): [number, number] + /** World square to WGS84 degrees; null where the CRS has no answer. */ + toLngLat(x: number, y: number): [number, number] | null +} + +const mappings = new WeakMap() + +/** + * How far, in world units, a round trip may land from where it started for + * the inverse to count, MapLibre's own tolerance: loose enough for a + * converter that is only approximately invertible, tight enough to reject an + * answer from the other side of the world. + */ +const ROUND_TRIP_WORLD = 1e-3 + +function worldMappingOf(definition: MapProjection): WorldMapping { + let mapping = mappings.get(definition) + if (mapping) return mapping + const { + projection, + tileMatrixSet: { + origin: [originX, originY], + extentAtZoom0: extent, + }, + } = definition + const toWorld = (lng: number, lat: number): [number, number] => { + if (!Number.isFinite(lng) || !Number.isFinite(lat)) return [NaN, NaN] + try { + const [x, y] = projection.forward([lng, lat]) + return Number.isFinite(x) && Number.isFinite(y) + ? [(x - originX) / extent, (originY - y) / extent] + : [NaN, NaN] + } catch { + return [NaN, NaN] + } + } + // An inverse that forward does not take back to the position has no + // answer there, as MapLibre treats it. + const toLngLat = (x: number, y: number): [number, number] | null => { + try { + const lngLat = projection.inverse([ + originX + x * extent, + originY - y * extent, + ]) + if ( + !lngLat || + !Number.isFinite(lngLat[0]) || + !Number.isFinite(lngLat[1]) + ) { + return null + } + const [backX, backY] = toWorld(lngLat[0], lngLat[1]) + return Math.hypot(backX - x, backY - y) <= ROUND_TRIP_WORLD + ? [lngLat[0], lngLat[1]] + : null + } catch { + return null + } + } + mapping = { toWorld, toLngLat } + mappings.set(definition, mapping) + return mapping +} + +/** + * The world of the map's projection when it is one registered with + * `addProjection` and among the projections given, or null otherwise. The + * same projection gives the same object each frame. + */ +export function worldMappingFor( + map: MapLike | null, + projections: readonly MapProjection[] +): WorldMapping | null { + const name = customProjectionName(map?.getProjection?.()) + if (!name) return null + const definition = projections.find((p) => p.name === name) + return definition ? worldMappingOf(definition) : null +} + +/** Source CRS to WGS84 and back, as the layer's projection context holds it. */ +export interface LngLatTransformer { + forward(x: number, y: number): [number, number] + inverse(lng: number, lat: number): [number, number] +} + +/** A box in the world square. */ +export interface WorldBounds { + x0: number + y0: number + x1: number + y1: number +} + +/** Source CRS to the world square, through WGS84; NaN where either has no answer. */ +export function sourceToWorld( + mapping: WorldMapping, + to4326: LngLatTransformer +): (x: number, y: number) => [number, number] { + return (x, y) => { + const [lng, lat] = to4326.forward(x, y) + return mapping.toWorld(lng, lat) + } +} + +/** World square to source CRS, through WGS84; NaN where either has no answer. */ +export function worldToSource( + mapping: WorldMapping, + to4326: LngLatTransformer +): (x: number, y: number) => [number, number] { + return (x, y) => { + const lngLat = mapping.toLngLat(x, y) + if (!lngLat) return [NaN, NaN] + const source = to4326.inverse(lngLat[0], lngLat[1]) + return Number.isFinite(source[0]) && Number.isFinite(source[1]) + ? source + : [NaN, NaN] + } +} + +/** + * The world-square box around a source-CRS rectangle, from a grid of + * samples across it, so a pole inside the rectangle counts. Its edges are + * sampled twice as densely, and the box grows by half the farthest any + * in-between sample bends from its neighbors, twice what a smooth edge can + * bulge between samples, so an edge that curves out between them stays + * inside. Null when no sample has a place in the world. + */ +export function worldBoundsOf( + bounds: { xMin: number; xMax: number; yMin: number; yMax: number }, + mapping: WorldMapping, + to4326: LngLatTransformer, + samples = 8 +): WorldBounds | null { + const toWorld = sourceToWorld(mapping, to4326) + let x0 = Infinity + let y0 = Infinity + let x1 = -Infinity + let y1 = -Infinity + const record = (wx: number, wy: number) => { + if (!Number.isFinite(wx) || !Number.isFinite(wy)) return + x0 = Math.min(x0, wx) + y0 = Math.min(y0, wy) + x1 = Math.max(x1, wx) + y1 = Math.max(y1, wy) + } + const at = (u: number, v: number) => + toWorld( + bounds.xMin + u * (bounds.xMax - bounds.xMin), + bounds.yMin + v * (bounds.yMax - bounds.yMin) + ) + for (let j = 0; j <= samples; j++) { + for (let i = 0; i <= samples; i++) { + record(...at(i / samples, j / samples)) + } + } + let bend = 0 + const steps = 2 * samples + for (const edge of [ + (t: number) => at(t, 0), + (t: number) => at(t, 1), + (t: number) => at(0, t), + (t: number) => at(1, t), + ]) { + const points = Array.from({ length: steps + 1 }, (_, k) => edge(k / steps)) + for (let k = 1; k < steps; k += 2) { + const [ax, ay] = points[k - 1] + const [mx, my] = points[k] + const [bx, by] = points[k + 1] + record(mx, my) + const offset = Math.hypot(mx - (ax + bx) / 2, my - (ay + by) / 2) + if (Number.isFinite(offset)) bend = Math.max(bend, offset) + } + } + if (!Number.isFinite(x0)) return null + const pad = bend / 2 + return { x0: x0 - pad, y0: y0 - pad, x1: x1 + pad, y1: y1 + pad } +} + +/** + * The share of the world square's width a source-CRS rectangle spans inside + * it, from a grid of samples across it that land in the square. Points past + * the square, such as a pole the projection sends far away, do not count. + * Undefined when the samples in the square span no width. + */ +export function worldFractionOf( + bounds: { xMin: number; xMax: number; yMin: number; yMax: number }, + mapping: WorldMapping, + to4326: LngLatTransformer, + samples = 32 +): number | undefined { + const toWorld = sourceToWorld(mapping, to4326) + let x0 = Infinity + let x1 = -Infinity + for (let j = 0; j <= samples; j++) { + const y = bounds.yMin + (j / samples) * (bounds.yMax - bounds.yMin) + for (let i = 0; i <= samples; i++) { + const x = bounds.xMin + (i / samples) * (bounds.xMax - bounds.xMin) + const [wx, wy] = toWorld(x, y) + if (!(wx >= 0 && wx <= 1 && wy >= 0 && wy <= 1)) continue + x0 = Math.min(x0, wx) + x1 = Math.max(x1, wx) + } + } + return x1 > x0 ? x1 - x0 : undefined +} diff --git a/src/zarr-layer.ts b/src/zarr-layer.ts index c598c8b1..8fe48738 100644 --- a/src/zarr-layer.ts +++ b/src/zarr-layer.ts @@ -39,8 +39,14 @@ import { computeWorldOffsets, resolveProjectionParams, isGlobeProjection as checkGlobeProjection, + customProjectionName, } from './map-utils' import { MAPBOX_IDENTITY_MATRIX } from './mapbox-utils' +import { + worldMappingFor, + type MapProjection, + type WorldMapping, +} from './world-mapping' import { mapboxEcefClipMatrix, parseMapboxRenderParameters, @@ -216,8 +222,10 @@ export class ZarrLayer { private onAuthError: OnAuthError | undefined private customStore: Readable | undefined private renderPoles: boolean + private mapProjections: readonly MapProjection[] private lastIsGlobe: boolean | null = null - private lastZoomInRange: boolean | null = null + private lastShouldDraw: boolean | null = null + private warnedProjection: string | null = null private usingDirectMapboxGlobePath: boolean = false private renderPathSwitchPending: boolean = false private mapboxRenderParametersSeen: boolean = false @@ -368,6 +376,7 @@ export class ZarrLayer { onAuthError, store, renderPoles = false, + mapProjections = [], }: ZarrLayerOptions) { if (!id) { throw new Error('[ZarrLayer] id is required') @@ -443,6 +452,7 @@ export class ZarrLayer { this.onAuthError = onAuthError this.customStore = store this.renderPoles = renderPoles + this.mapProjections = mapProjections } /** @@ -679,7 +689,8 @@ export class ZarrLayer { this.lastIsGlobe = isGlobe this.regionRenderer?.onProjectionChange(isGlobe) - if (this.isZoomInRange()) this.regionRenderer?.update(this.map, this.gl!) + this.syncWorldMapping() + if (this.shouldDraw()) this.regionRenderer?.update(this.map, this.gl!) else this.regionRenderer?.followZoom(this.map, this.gl!) } catch (err) { this.initError = err instanceof Error ? err : new Error(String(err)) @@ -731,7 +742,8 @@ export class ZarrLayer { // per-frame update picks up again once it zooms back in. Until then the // level still follows the zoom, for queries and `ready`. if (this.map && this.gl) { - if (this.isZoomInRange()) this.regionRenderer.update(this.map, this.gl) + this.syncWorldMapping() + if (this.shouldDraw()) this.regionRenderer.update(this.map, this.gl) else this.regionRenderer.followZoom(this.map, this.gl) } } @@ -822,6 +834,36 @@ export class ZarrLayer { return zoom >= this.minZoom && zoom <= this.maxZoom } + /** + * Whether the layer fetches and draws this frame: inside its zoom range and + * in a projection it can place data in. A projection registered with + * MapLibre's `addProjection` has to be among `mapProjections`. + */ + private shouldDraw(): boolean { + const customProjection = customProjectionName(this.map?.getProjection?.()) + const unplaceable = + customProjection && !this.worldMapping() ? customProjection : null + if (unplaceable && unplaceable !== this.warnedProjection) { + console.warn( + `[ZarrLayer] The map is in the ${unplaceable} projection, which is ` + + `not among the layer's mapProjections, so the layer is hidden ` + + `while it is active.` + ) + } + this.warnedProjection = unplaceable + return !unplaceable && this.isZoomInRange() + } + + /** The world of the map's projection, when it is a registered one. */ + private worldMapping(): WorldMapping | null { + return worldMappingFor(this.map, this.mapProjections) + } + + /** Place regions in the world of the map's current projection. */ + private syncWorldMapping(): void { + this.regionRenderer?.setWorldMapping(this.worldMapping()) + } + prerender( _gl: WebGL2RenderingContext | WebGLRenderingContext, _params: unknown, @@ -840,14 +882,15 @@ export class ZarrLayer { // Draped terrain tiles are cached by the map and reused across zoom // changes while it moves, so crossing the zoom range has to force a // redraw or the tiles keep (or lack) the layer until movement stops. - const inRange = this.isZoomInRange() - if (this.lastZoomInRange !== null && this.lastZoomInRange !== inRange) { + const draws = this.shouldDraw() + if (this.lastShouldDraw !== null && this.lastShouldDraw !== draws) { this.invalidate() } - this.lastZoomInRange = inRange - if (!inRange) return + this.lastShouldDraw = draws + if (!draws) return this.syncProjectionState() + this.syncWorldMapping() const cover = this.coveringTiles() // Tiles the draped hooks asked for are only being drawn while the map // drapes the layer: MapLibre with terrain, Mapbox with terrain or on its @@ -939,7 +982,7 @@ export class ZarrLayer { } this.noteMapboxRenderParameters(renderParameters) - if (!this.isZoomInRange()) { + if (!this.shouldDraw()) { return } @@ -977,7 +1020,11 @@ export class ZarrLayer { : undefined const isGlobe = this.isGlobeProjection() - const worldOffsets = computeWorldOffsets(this.map, isGlobe) + // A registered projection's world is drawn once, without copies. + const worldOffsets = computeWorldOffsets( + this.map, + isGlobe || !!this.worldMapping() + ) const colormapTexture = this.colormap.ensureTexture(this.gl) // The direct globe path draws straight onto the sphere, so it needs the // per-frame globe-world-to-clip matrix; without one this frame it falls @@ -1070,7 +1117,7 @@ export class ZarrLayer { return } - if (!this.isZoomInRange()) { + if (!this.shouldDraw()) { return } diff --git a/src/zarr-layer.zoom-range.test.ts b/src/zarr-layer.zoom-range.test.ts index cb3b23a0..675ad4ae 100644 --- a/src/zarr-layer.zoom-range.test.ts +++ b/src/zarr-layer.zoom-range.test.ts @@ -3,6 +3,7 @@ import { ZarrLayer } from './zarr-layer' import { buildMemoryZarrStore } from './__fixtures__/memory-zarr' import { createRecordingGl } from './__fixtures__/fake-gl' import type { MapLike } from './types' +import type { MapProjection } from './world-mapping' /** * A layer outside its zoom range must not fetch. The per-frame path checks @@ -90,10 +91,10 @@ function countingStore(times?: number, pyramid = false) { } } -function fakeMap(zoom: number) { - const state = { zoom } +function fakeMap(zoom: number, projection = 'mercator') { + const state = { zoom, projection } const map = { - getProjection: () => ({ type: 'mercator' }), + getProjection: () => ({ type: state.projection }), getTerrain: () => null, getZoom: () => state.zoom, getBounds: () => ({ @@ -119,10 +120,12 @@ async function addLayer( times?: number pyramid?: boolean selector?: ConstructorParameters[0]['selector'] + projection?: string + mapProjections?: MapProjection[] } = {} ) { const { store, chunkReads } = countingStore(opts.times, opts.pyramid) - const { map, state } = fakeMap(zoom) + const { map, state } = fakeMap(zoom, opts.projection) const gl = createRecordingGl() let ready: () => void @@ -142,6 +145,7 @@ async function addLayer( latIsAscending: false, minzoom, selector: opts.selector, + mapProjections: opts.mapProjections, onLoadingStateChange: (loadingState) => { if (!loadingState.metadata) ready() }, @@ -217,3 +221,58 @@ describe('zoom range', () => { expect(chunkReads.every((key) => key.includes('/c/1/'))).toBe(true) }) }) + +describe('custom projection', () => { + it('fetches in a registered projection the layer was given', async () => { + const { chunkReads } = await addLayer(6, 0, { + projection: 'custom:plate', + mapProjections: [ + { + name: 'custom:plate', + projection: { + forward: ([lng, lat]) => [lng, lat], + inverse: ([x, y]) => [x, y], + }, + tileMatrixSet: { origin: [-180, 90], extentAtZoom0: 360 }, + }, + ], + }) + expect(chunkReads.length).toBeGreaterThan(0) + }) + + it('fetches nothing and warns once in a registered projection it was not given', async () => { + const warn = vi.spyOn(console, 'warn').mockImplementation(() => {}) + try { + const { layer, gl, chunkReads } = await addLayer(6, 0, { + projection: 'EPSG:3031', + }) + layer.prerender(gl, {}) + layer.prerender(gl, {}) + await new Promise((resolve) => setTimeout(resolve, 20)) + + expect(chunkReads).toEqual([]) + const warnings = warn.mock.calls.filter(([message]) => + String(message).includes('EPSG:3031') + ) + expect(warnings).toHaveLength(1) + } finally { + warn.mockRestore() + } + }) + + it('fetches once the map returns to a built-in projection from one it was not given', async () => { + const warn = vi.spyOn(console, 'warn').mockImplementation(() => {}) + try { + const { layer, gl, state, chunkReads } = await addLayer(6, 0, { + projection: 'EPSG:3031', + }) + expect(chunkReads).toEqual([]) + + state.projection = 'mercator' + layer.prerender(gl, {}) + await until(() => chunkReads.length > 0) + } finally { + warn.mockRestore() + } + }) +})