@carbonplan/zarr-layer

Custom layer for rendering Zarr datasets in MapLibre or Mapbox GL, inspired (and borrowing significant code and concepts from) zarr-gl, zarr-cesium, @carbonplan/maps, and deck-gl-raster. Uses CustomLayerInterface to render data directly to the map and supports rendering to globe and mercator projections for both MapLibre and Mapbox. Input data are reprojected on the fly.

demo

See the demo for a quick tour of capabilities. Code for the demo is in the /demo folder.

data requirements

Supports v2 and v3 zarr stores via zarrita. Arbitrary CRS support via proj4 reprojection.

Self-describing stores

A store that carries the zarr proj and spatial conventions needs no configuration: the layer reads its CRS, extent, orientation and axis names straight out of the metadata, with no coordinate-array reads at all.

| Attribute | What it settles | | ----------------------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------- | | proj:code | The CRS by identifier: EPSG:4326, EPSG:3857, OGC:CRS84, the WGS84 UTM zones (EPSG:326xx/EPSG:327xx), or any code registered with proj4.defs. | | proj:wkt2 / proj:projjson | The CRS in full, for codes proj4 doesn't ship. No lookup table needed. | | spatial:transform | The extent and which edge row 0 sits on. | | spatial:bbox | The extent, winning over the transform's when both are declared. | | spatial:registration | Whether the declared coordinates fall on cell edges (pixel, the default) or cell centers (node). | | spatial:dimensions | Which array dimensions are y and x, for axes not named something recognizable. | | spatial:shape | Each pyramid level's size, on multiscales.layout entries, sparing an array open per level. |

Declared attributes are authoritative. If a store publishes them, they are used as-is rather than checked against its coordinate arrays, so a store whose attributes disagree with its own data will render wrong rather than be quietly corrected.

Constructor options always win over what the store declares: crs, proj4, bounds, latIsAscending and spatialDimensions each override the corresponding attribute.

For a store declaring a proj:code proj4 doesn't know and no proj:wkt2 or proj:projjson, the layer warns and leaves the CRS unresolved. Pass the definition as the proj4 prop, or register the code with proj4 to have the store's own declaration resolve (see registering a CRS).

Where these attributes are absent the layer falls back to what it always did: reading the coordinate arrays for extent and orientation, matching dimension names against a list of common aliases, and inferring the CRS from the magnitude of the bounds. A CF grid_mapping variable is read for crs_wkt when no proj: attribute is present.

Multiscales

High resolution datasets require multiscales. Chunks are loaded based on viewport intersection, and the level is chosen to match the screen resolution. Supports the zarr multiscales convention and legacy ndpyramid outputs, and tries to interpret other multiscale formats. See topozarr for a look at how to create these datasets.

globe rendering and polar coverage

Web Mercator rendering clips near ±85° latitude, leaving visible "pole holes" on globe projections. For EPSG:4326 and proj4 datasets:

MapLibre — Full polar coverage is always enabled via a direct ECEF rendering path. No configuration needed.

Mapbox — Set renderPoles: true to enable an experimental direct ECEF path that bypasses tile draping. This relies on Mapbox internal APIs and may break across Mapbox GL JS versions. During Mapbox's globe-to-mercator zoom morph the layer automatically falls back to the standard draped path (with pole holes visible). Incompatible with draping the zarr layer over Mapbox terrain — when terrain is enabled the layer always uses the draped tile path.

new ZarrLayer({
  // ...
  renderPoles: true, // Mapbox only — MapLibre always renders to the poles
})

install

npm install @carbonplan/zarr-layer

build locally

npm install
npm run build

usage

import maplibregl from 'maplibre-gl' // or mapbox
import { ZarrLayer } from '@carbonplan/zarr-layer'

const map = new maplibregl.Map({container: 'map'})
const layer = new ZarrLayer({
  id: 'zarr-layer',
  source: 'https://example.com/my.zarr',
  variable: 'temperature',
  clim: [270, 310],
  colormap: ['#000000', '#ffffff', ...],
  selector: { month: 1 },
})
map.on('load', () => {
  map.addLayer(layer)
  // optionally add before id to slot data into map layer stack.
  // map.addLayer(layer, 'beforeID')
})

options

Required: | Option | Type | Description | |--------|------|-------------| | id | string | Unique layer identifier | | source | string | Zarr store URL (required unless store is provided) | | variable | string | Variable name to render | | colormap | array | Array of hex strings or [r,g,b] values | | clim | [min, max] | Color scale limits |

Optional: | Option | Type | Default | Description | |--------|------|---------|-------------| | store | Readable | - | Custom zarrita-compatible store (e.g., IcechunkStore). When provided, source becomes optional. | | selector | object | {} | Dimension selector (unspecified dims default to index 0) | | opacity | number | 1 | Layer opacity (0-1) | | zarrVersion | 2 | 3 | auto | Zarr format version (tries v3 first, falls back to v2) | | minzoom | number | 0 | Minimum zoom level for rendering | | maxzoom | number | Infinity | Maximum zoom level for rendering | | fillValue | number | auto | No-data value (from metadata if not set) | | spatialDimensions | object | auto | Custom { lat, lon } dim names | | crs | string | auto | CRS identifier. Not needed for EPSG:4326/EPSG:3857 data (detected automatically). Codes proj4 defines (the WGS84 UTM zones, among others) or that were registered with proj4.defs work without a proj4 string. | | proj4 | string | - | Proj4 definition string for CRS reprojection (bounds recommended, else derived from coordinates) | | bounds | array | auto | [xMin, yMin, xMax, yMax] in source CRS units (degrees for EPSG:4326, meters for EPSG:3857). These are interpreted as edge bounds (not center-to-center) | | latIsAscending | boolean | auto | Latitude orientation | | renderingMode | '2d' | '3d' | '3d' | Custom layer rendering mode | | customFrag | string | - | Custom fragment shader | | uniforms | object | - | Shader uniform values (requires customFrag) | | onLoadingStateChange | function | - | Loading state callback | | transformRequest | function | - | Transform request URLs and add headers/credentials (see authentication) | | onAuthError | function | - | Called with the HTTP status when a signed request fails with expired credentials (see authentication) | | renderPoles | boolean | false | Enable polar coverage in Mapbox globe for EPSG:4326/proj4 datasets (see globe rendering). No effect on EPSG:3857 data. MapLibre always renders to the poles. |

methods

layer.setOpacity(0.8)
layer.setClim([0, 100])
layer.setColormap(['#000', '#fff'])
layer.setSelector({ time: 5 })
layer.setVariable('precipitation') // async - reloads metadata
layer.setUniforms({ u_weight: 1.5 }) // no-op unless layer has customFrag

throttling rapid selector changes

setSelector is synchronous and kicks off a fetch on every call. For UIs that fire rapid updates (e.g. dragging a time slider), debounce the value on your side before handing it to the layer so you aren't firing one abort-cancelled fetch per pointer event.

selectors

Selectors specify which slice of your multidimensional data to render. Dimensions not specified default to index 0.

Basic syntax:

// Simple value - matches exact value in coordinate array
{ time: 5 }
{ time: '2024-01-15' }

// Explicit index - uses array index directly (no coordinate lookup)
{ time: { selected: 5, type: 'index' } }

// Explicit value - same as simple syntax, matches exact value
{ time: { selected: 5, type: 'value' } }

Multi-band selection (for custom shaders):

// String values use the string directly as the shader variable name
{ band: ['tavg', 'prec'] }
// exposes as: tavg, prec

// Numeric values are prefixed with the dimension key (required for valid GLSL identifiers)
{ month: [1, 2, 3] }
// exposes as: month_1, month_2, month_3

// Index selectors use dimension-prefixed shader variable names
{ band: { selected: [0, 1, 2], type: 'index' } }
// exposes as: band_0, band_1, band_2; update shader references accordingly

// Mix with other dimensions
{ band: ['red', 'green', 'blue'], time: 0 }

Query-specific selectors:

// Array of values for time series queries
const result = await layer.queryData(
  { type: 'Point', coordinates: [lng, lat] },
  { time: [0, 1, 2, 3, 4] } // returns data for all 5 time steps
)

Type options:

| Type | Behavior | | ------------------- | ------------------------------------------------------------- | | 'value' (default) | Matches exact value in coordinate array (throws if not found) | | 'index' | Uses value directly as array index |

custom shaders and uniforms

Custom fragment shaders let you do math on your data to change how it's displayed. This can be useful for things like log scales, combining bands, or aggregating data over a time window. Bands can span separate chunks — each band is fetched in parallel and combined for rendering. You can pass in uniforms to allow user interaction to influence the custom shader code.

Band names are automatically sanitized to valid GLSL identifiers: any characters that aren't letters, digits, or underscores are replaced with underscores, and names starting with a digit are prefixed with an underscore. For example, s2med_harvest:B02 becomes s2med_harvest_B02 and 123band becomes _123band.

The layer blends with a premultiplied alpha blend function, so your customFrag must output premultiplied color. Multiply RGB by your final alpha (e.g. fragColor = vec4(c.rgb * opacity, opacity)). Emitting straight vec4(c.rgb, opacity) renders correctly only at full opacity.

When a customFrag is supplied, it owns discarding. Missing/fill pixels are surfaced as NaN instead of being dropped automatically, so you can aggregate over partial coverage (e.g. a mean over bands with differing coverage) rather than just their intersection. Note that NaN propagates through arithmetic (even NaN * 0.0 is NaN), so guard values before multiplying or accumulating, e.g. isnan(x) ? 0.0 : x. Add your own discard for pixels you want to drop:

new ZarrLayer({
  // ...
  customFrag: `
    uniform float u_weight;
    if (isnan(band_a)) {
      discard;
    }
    float val = band_a * u_weight;
    float norm = (val - clim.x) / (clim.y - clim.x);
    vec4 c = texture(colormap, vec2(clamp(norm, 0.0, 1.0), 0.5));
    fragColor = vec4(c.rgb * opacity, opacity);
  `,
  uniforms: { u_weight: 1.0 },
})

NDVI example

Here's an example of computing NDVI (Normalized Difference Vegetation Index) using custom shaders:

new ZarrLayer({
  source: 'https://example.com/sentinel2.zarr',
  variable: 'data',
  colormap: [
    /* gradient */
  ],
  selector: { band: ['B08', 'B04'], time: 0 },
  clim: [-1, 1],
  customFrag: `
    if (isnan(B08) || isnan(B04)) {
      discard;
    }
    float ndvi = (B08 - B04) / (B08 + B04);
    float norm = (ndvi - clim.x) / (clim.y - clim.x);
    vec4 c = texture(colormap, vec2(clamp(norm, 0.0, 1.0), 0.5));
    fragColor = vec4(c.rgb * opacity, opacity);
  `,
})

custom projections

Datasets in EPSG:4326 or EPSG:3857 need no CRS configuration. For anything else (e.g., Lambert Conformal Conic, UTM), set crs to the code. It resolves with no further configuration if proj4 defines that code or you registered it yourself; otherwise pass a proj4 definition string alongside it, or the renderer will warn and fall back to inferred CRS. Specifying bounds in source CRS units is recommended for performance (otherwise derived from coordinate arrays).

None of this is needed for a store carrying the proj and spatial conventions (see self-describing stores).

new ZarrLayer({
  // ...
  spatialDimensions: {
    lat: 'projection_y_coordinate',
    lon: 'projection_x_coordinate',
  },
  proj4:
    '+proj=lcc +lat_1=38.5 +lat_2=38.5 +lat_0=38.5 +lon_0=-97.5 +x_0=0 +y_0=0 +R=6371229 +units=m +no_defs',
  bounds: [-2697520, -1587306, 2697480, 1586694], // recommended: edge bounds [xMin, yMin, xMax, yMax] in source CRS units
})

The data will be reprojected to Web Mercator for display using GPU-accelerated mesh reprojection powered by @developmentseed/raster-reproject. Find proj4 strings at epsg.io or in your dataset's metadata.

registering a CRS

proj4 defines only a small set of codes out of the box: EPSG:4326, EPSG:4269, EPSG:3857 (with its aliases), EPSG:5041, EPSG:5042, and the 120 WGS84 UTM zones (EPSG:32601-EPSG:32660 north, EPSG:32701-EPSG:32760 south). Everything else needs a definition from you, including national grids like EPSG:27700 and EPSG:2154, and the UTM zones on other datums, which look built-in but aren't: NAD83 (EPSG:269xx) and ETRS89 (EPSG:258xx) are separate codes from their WGS84 counterparts.

Pass it as the proj4 prop, which takes precedence over crs and over anything the store declares. Or register it with proj4 before creating the layer, which lets a store's own proj:code resolve with no per-layer configuration:

import proj4 from 'proj4'

proj4.defs('EPSG:25833', '+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs')

Registering pays off only when you don't know which store is in which CRS ahead of time; otherwise the prop is simpler. proj4.defs also accepts WKT2 and PROJJSON. Note that zarr-layer imports proj4 rather than bundling its own copy, so your proj4 import is the registry it reads, as long as your bundler resolves both to one copy.

A store declaring proj:wkt2 or proj:projjson carries its own definition and needs none of this.

queries

Supports Point, Polygon, and MultiPolygon geometries in geojson format. You can optionally pass in a custom selector to override the visualization selector.

// Point query
const result = await layer.queryData(
  { type: 'Point', coordinates: [lng, lat] },
  // optional selector override (useful for e.g. time series creation)
  { time: [0, 1, 2] }
)

// Polygon query
const result = await layer.queryData({
  type: 'Polygon',
  coordinates: [[...]],
})

// Returns:
// {
//   [variable]: number[],
//   dimensions: ['<store-y-axis>', '<store-x-axis>'],
//   coordinates: { '<store-y-axis>': number[], '<store-x-axis>': number[] }
// }

Spatial query results are returned in the dataset's source CRS, under the store's own spatial axis names. An EPSG:3857 dataset with y/x axes returns Web Mercator meters under y/x; an EPSG:4326 dataset with latitude/longitude axes returns degrees under latitude/longitude; a custom-proj4 dataset returns its source-CRS values under whatever the store calls them. Input geometries are still supplied as GeoJSON lon/lat regardless of the source CRS.

You can pass a third options argument to control query behavior:

const result = await layer.queryData(geometry, selector, {
  signal: abortController.signal, // cancel in-flight query
  includeSpatialCoordinates: false, // omit per-pixel coordinates for slimmer results
  level: 'finest', // read the highest-resolution level instead of the drawn one
})

Note: Query results match rendered values (scale_factor/add_offset applied, fillValue/NaN filtered).

query resolution

By default a query reads the level the map is currently drawing, so results agree with what the user sees and zooming out coarsens them. Pass level: 'finest' to always read the highest-resolution level in the store, which is what you want when the answer shouldn't depend on the camera — sampling point features, for instance.

'finest' reads a level the renderer may not hold, so it fetches cold instead of reusing chunks the render path already cached. A point costs about one chunk either way; a polygon covers quadratically more pixels at a finer level, so on a deep pyramid at low zoom it can read many more. It doesn't disturb rendering: the query reads its own level and leaves the drawn one alone. No effect on single-level stores.

query readiness

queryData waits for metadata and a committed resolution level, so it can be called immediately after map.addLayer(layer) with no render pass in between and no polling. Readiness failures — initialization failure, failure to load a level, removal from the map, or querying before the layer was added — reject with a ZarrLayerNotReadyError rather than returning empty. Initialization and level-load failures carry the underlying error on .cause.

Failed reads reject as well, so an empty result means the geometry found no data.

layer.ready exposes the same wait as a promise, for uses other than queries:

map.addLayer(layer)
await layer.ready // metadata loaded and a resolution level committed

Not the same signal as onLoadingStateChange, which is a spinner: it flips as chunks load and reports nothing about the level commit, so loading: false can be emitted while the layer still has no level. Its error field reports initialization and level-load failures and returns to null after the load recovers or the selector changes.

authentication

Use transformRequest to add headers or credentials to requests. The function receives the fully resolved URL for each request, enabling per-path authentication like presigned S3 URLs. Supports any fetch options.

// Static auth (same headers for all requests)
transformRequest: (url) => ({
  url,
  headers: { Authorization: `Bearer ${token}` },
})

// Presigned URLs (path-specific signatures)
transformRequest: async (url) => ({
  url: await getPresignedUrl(url),
})

expired credentials

Credentials signed into a request expire mid-session. Those failures are otherwise invisible, because the layer treats a rejected read as a missing chunk and renders a hole. onAuthError surfaces the status so you can refresh and re-add the layer:

onAuthError: async (status) => {
  await refreshCredentials()
  map.removeLayer(layer.id)
  map.addLayer(newLayerWithFreshCredentials())
}

It fires on 400 and 401 only, and at most once per store, so a burst of concurrent chunk failures triggers a single refresh. 400 is included because expired temporary AWS credentials return ExpiredToken/InvalidToken as a 400, often with no readable body on HEAD probes or CORS-gated responses. 403 is excluded: S3 and CloudFront return it for genuinely absent chunks, which sparse pyramids hit routinely. Without a handler, a 400 propagates as an error rather than being read as a missing chunk.

custom stores

For advanced use cases like Icechunk, you can pass a custom zarrita-compatible store directly. When using a custom store, source becomes optional:

import { IcechunkStore } from 'icechunk-js'

const store = await IcechunkStore.open(...)

new ZarrLayer({
  id: 'icechunk-layer',
  store,
  variable: 'temperature',
  colormap: [...],
  clim: [0, 100],
})

The store must implement the zarrita Readable interface with at minimum a get(key: string) method.

codecs

The library uses zarrita for Zarr data access. zarrita includes built-in codecs for bytes, zlib, gzip, blosc, lz4, zstd, transpose, crc32c, and bitround. You can add more if needed!

adding custom codecs

Virtualized NetCDF data may use numcodecs.*-prefixed codec names (e.g., numcodecs.zlib, numcodecs.shuffle) that zarrita doesn't recognize by default. Use codecRegistry to register them before creating layers:

import { codecRegistry } from '@carbonplan/zarr-layer'

// Alias numcodecs.zlib → zarrita's built-in zlib
const zlibFactory = codecRegistry.get('zlib')
if (zlibFactory) codecRegistry.set('numcodecs.zlib', zlibFactory)

// numcodecs.shuffle — byte un-shuffle by element size
codecRegistry.set('numcodecs.shuffle', async () => ({
  fromConfig(config: { elementsize?: number }) {
    const elementsize = config?.elementsize ?? 1
    return {
      kind: 'bytes_to_bytes',
      decode(bytes: Uint8Array): Uint8Array {
        if (elementsize <= 1) return bytes
        const n = bytes.length
        const count = Math.floor(n / elementsize)
        const out = new Uint8Array(n)
        for (let i = 0; i < count; i++) {
          for (let j = 0; j < elementsize; j++) {
            out[i * elementsize + j] = bytes[j * count + i]
          }
        }
        for (let i = count * elementsize; i < n; i++) {
          out[i] = bytes[i]
        }
        return out
      },
    }
  },
}))

thanks

This experiment is only possible following in the footsteps of other work in this space. zarr-gl showed that custom layers are a viable rendering option and zarr-cesium showed how flexible web rendering can be. We borrow code and concepts from both. This library also leans on our prior work on @carbonplan/maps for many of its patterns. Custom projection support uses @developmentseed/raster-reproject for adaptive mesh generation. LLMs of several makes aided in the coding and debugging of this library.

license

All the code in this repository is MIT-licensed, but we request that you please provide attribution if reusing any of our digital content (graphics, logo, articles, etc.).

about us

CarbonPlan is a nonprofit organization that uses data and science for climate action. We aim to improve the transparency and scientific integrity of climate solutions with open data and tools. Find out more at carbonplan.org or get in touch by opening an issue or sending us an email.