Band math over time, step by step.

01The idea

A calculator in a GIS combines rasters cell by cell and writes a file before you see anything. Janus computes as it reads: an expression becomes a small program, compiled once and evaluated over blocks of 4096 pixels wherever the series is read, the map's overview, the tiles you zoom into, a pixel's series under the cursor, a region, a Zeit tool. Nothing is written until you ask for it.

An index

On the Display panel, the series shown as an index of its bands: EVI, SAVI, MNDWI, NBR, any of the 275 of the Awesome Spectral Indices the series has the bands for, or a formula of your own.

The calculator

Tools → Raster calculator...: an expression over every open layer, each by a short name, a series with an image, two series matched by date, a series reduced over time, rules. The result is a new layer.

02The data

The series of the STAC guide: Planetary Computer's Sentinel-2 L2A over a square of central Rondônia, Brazil, 8.8 km a side, the 107 dates from June 2019 to September 2024 whose scenes had 5% of clouds at most (most of them in the dry seasons), the clouds left masked by the scene classification (SCL). Here it was downloaded as GeoTIFFs: blue, green, red, NIR, the two SWIR and SCL, 877 × 886 cells of 10 m. The calculator works the same on the series read in place from the catalog; only the reading is slower.

Two pins follow the guide: pin 1 in a forest in the north-east (its NDVI around 0.8 every year), pin 2 in a pasture in the west (green in the wet months, lower every dry season: 0.18 at its lowest, in September 2024).

03An index on the series

The Display panel has, under the bands, Index or expression. Its list holds the indices of the Awesome Spectral Indices (Montero et al. 2023) the series can compute, searched by name or domain (vegetation, water, burn, urban, soil, snow). Each index names its bands by role: N the NIR, R the red, S1 the first SWIR, RE1 the first red edge. Janus finds the series' band for each role by its wavelength or its name (B08, NIR, nir08, SR_B5), so an index of Sentinel-2 works on Landsat too. Its constants take the catalogue's defaults (L = 0.5 for SAVI), its wavelengths (lambdaN) those of the series' bands.

The formula can be edited, or typed: 2.5 * (N - R) / (N + 6*R - 7.5*B + 1), nd(nir, swir1), (B08 - B04) / (B08 + B04). It is checked as you type, and Apply shows the series as the expression: the map at every date and in every mode, the chart, the pins, the region, the exports, the Zeit tools and Copy Python, which writes the same formula in numpy. From the command line: jn SERIES --expr EVI, or --expr "(N - S1) / (N + S1)".

04The calculator

Tools → Raster calculator... (or Open in calculator from the Display panel) opens a window with an editor on the left and, on the right, four lists: the open Layers with their short names and bands, the Functions, the Indices and the History.

Short names

Each open layer gets one from its name (s2, dem, l8); a band is s2.nir, s2.B08, s2.band4 or s2[4], a role s2.N. Bare names (N, nir, B04) are the bands of the active layer, or of the layer chosen under Options.

As you type

The text is coloured by what each name is (a layer, a band, a function, a variable of your own); a list completes names (Tab or Enter, Ctrl+Space to call it); the pointer over a name says what it is (a layer's dates and grid, a band's wavelength, a function's arguments); an error is underlined where it is, and clicking its message puts the cursor there.

Before running

Under the editor, what the result will be: a series of 107 dates or an image, its grid, its CRS, the layers read.

Run

Ctrl+Enter or Run. The text goes to the history; a star keeps it among the favourites.

05The language

One statement a line (or apart by ;), # for comments. A line name = expression names a value for the lines after it; the last line is the result.

No data is not a number that spreads: arithmetic on it gives no data, a reduction skips it, a comparison with it is no data, and coalesce or a rule x == nan = 0 fills it.

06Over time

A series has dates. An expression of a series is a series: (N - R) / (N + R) is a new value at each of the 107 dates. A reduction with one argument turns it into an image, over every date with a value:

# how many of the 107 dates had little vegetation
ndvi = (N - R) / (N + R)
count(ndvi < 0.4)

With more than one argument the same functions work pixel by pixel at each date: max(B04, B03, B02) is the brightest visible band of each date. A series and its reduction mix: ndvi - mean(ndvi) is each date's departure from the pixel's mean, a series again; (ndvi - mean(ndvi)) / std(ndvi) its z-score.

Time is an axis of the language, counted in the series' real dates (an irregular series with clouds stays irregular, a masked date does not exist):

The date of a value

argmax(ndvi) is the date of the pixel's largest NDVI (the peak of the season) as a decimal year, argmin of the smallest; first(ndvi < 0.3 and year >= 2020) and last(...) the first and last date a condition held. A result in years opens with a legend of years.

A slice of the calendar

median(ndvi[2019-06 .. 2019-09]) reduces June to September 2019 alone; max(ndvi[2019 .. 2021]) three years; mean(ndvi[11 .. 2]) November to February of every year.

By year, by month

max(ndvi) by year is a new series with one date per year, mean(ndvi) by month one per month that has a date; they mix with each other (max(ndvi) by year - min(ndvi) by year), not with the dates they came from.

A window in days

mean(ndvi, 90d) is, at each date, the mean of the dates of the 90 days up to it (pandas' rolling("90D")); w counts weeks, y years.

The value before

ndvi[-1] (or lag(ndvi, 1)) is the pixel's previous observation with a value, so ndvi - ndvi[-1] is the change since; lag(ndvi, 1y) the value of the nearest date a year before (within half the usual step between dates).

The anomaly

anomaly(ndvi, base=2015..2019, by=month) is each date's z-score against the mean and the deviation of its month in the base years; by=doy, window=15d against the days of the year within 15 days.

Events

count(ndvi < 0.3) the dates below; longest(ndvi < 0.3) the longest it stayed below, in days; persist(ndvi < 0.3, 3) a series, 1 where the last three observations were all below.

Tick Live and the result is on the map as it is typed, at the overview's size, as a layer marked preview: a reduction is drawn as its dates are read (max so far: 23 of 107 dates), a series at the date shown. Tick Chart and under the editor are its values at the cursor (or the first pin): the result's series, or for an image its value as a level (a date as a time) with the values its reductions reduce; select a part of the text and its series is drawn too, so first(ndvi < 0.3) is seen falling on the date the NDVI crossed 0.3.

07Rules

A line condition = value changes the pixels where the condition holds and keeps the others. The result starts as the condition's first operand, or as what result = ... made it, and each rule works on what the rules before left:

# the median NDVI of the five years, in three classes
m = median(nd(N, R))
result = m
m < 0.45 = 1
m >= 0.45 && m < 0.7 = 2
m >= 0.7 = 3

Rules are short for where: img > 0 = 1 is where(img > 0, 1, img). So img == nan = 0 fills the gaps with 0, result == 0 = nan makes 0 no data, and img.band1 > 1 && img.band2 >= 5 = 0 zeroes a stack's pixels where two of its bands agree.

08Several layers

The expression may use any open layer, in any grid or CRS:

A series and an image

The image is used at every date: s2.nir / dem, ndvi / max(ndvi).

Two series

Matched by date (Options → Match dates): the same day, the nearest within some days (Sentinel-2 with Landsat), the same month or the same year. A date with no match is no data there.

Another grid

A layer of other cells (Landsat's 30 m with Sentinel-2's 10 m) is resampled onto the result's grid, a layer in another CRS reprojected; the result takes the grid of the layer chosen under Grid of, else that of the series whose dates it has.

A stack

One file of several bands that are not dates (a DEM with its slope, a classification with its probabilities) is an image whose bands are named: terrain.slope, terrain.band2.

09The result

The result is a layer like any other: the map, the chart, pins, regions, the Zeit tools, exports. It is written in a folder of its own under %LOCALAPPDATA%\Janus\results\calc, as one of two things:

A layer of VRTs

The default: one small VRT per date whose band is computed by Janus's pixel function from the layers' files when it is read. It costs nothing to make; the map computes the overview and the tiles as it shows them.

GeoTIFFs

Write GeoTIFFs: the VRTs read once into cloud-optimized GeoTIFFs (DEFLATE, overviews), several dates at a time, as a task in Tools → Tasks. For a result to keep, to share or to open in another program.

The folder keeps the text and the layers it was computed from (calc.json) and the task's log. From the command line, jn SERIES --calc "max(nd(N, R))" --calc-run opens the series, runs the text and adds the result.

10Checked against the originals

The calculator was run against the Python libraries it stands in for, on the same files, comparing the values and the time each took (tools/bench_calc_vs_originals.py; the numbers in bench/calc-vs-originals.json). It was run on a 20-thread laptop held to three of its threads (Janus's reads and NumPy's alike), so the times move by about a tenth from one run to the next.

The indices, against spyndex 0.12

spyndex computes the same Awesome Spectral Indices catalogue with NumPy. All 184 indices that both can compute from six bands (two more need a second NIR band that spyndex does not default) were run over 2048 × 2048 random reflectances. Janus matched spyndex within 3 × 10⁻⁷ relative to its value (MSAVI the farthest, the rounding of single precision) and had no data on the same pixels. Together they took 3.7 s in Janus on one thread against 6.6 s in spyndex.

The reductions, against xarray 2025.12

A series of 48 monthly dates of 1024 × 1024 pixels with a fifth of every date clouded (NaN). xarray read the files with rasterio, Janus through the result's VRT as a GeoTIFF translation reads it; the times include the reading.

Over time, against xarray and pandas

An irregular series of 64 dates 16 to 20 days apart over three years, 512 × 512 pixels with a fifth of every date clouded. Janus's times include writing the result (the reductions read by several dates computed once into GeoTIFFs, as when the result is written out) and reading every date of it.

The dates in years differ by single precision (2021.5 is kept to about 10⁻⁴ of a year), the anomaly's z-scores likewise; where a month's base had no deviation xarray gives an infinity and Janus no data, and those pixels were left out. xarray holds the whole series in memory and subtracts it once; Janus writes and reads each of the 64 dates as files, so a series returned whole gains least, a reduction most.

Every result had no data on the same pixels as xarray's. A count is a sum of comparisons, cheap for both; its time is mostly reading the 48 files, which Janus does on its threads, each block of a file read once. xarray's quantile with NaNs sorts the values of every pixel one by one; Janus selects the rank in place.

On the map, an expression costs little over a band. Reading the overview of 24 dates of 2048 × 2048 pixels took 28 ms for a band, 66 ms for NDVI as a normalized difference, 81 ms for NDVI written as an expression, and 122 ms for EVI, whose three bands are read at every date. On the five years of Sentinel-2 of this guide, a count over its 107 dates like section 06's (count(nd(N, R) > 0.5)) is on the map in 0.8 s on three threads, 2.6 s when GDAL's VRT reads it; a pin's value at the 107 dates takes 0.1 s, and the median written as a cloud-optimized GeoTIFF 1.1 s (3.7 s). Janus reads its VRTs itself: each date's file straight, each block of a file by one thread, the values those of GDAL's VRT bit for bit; GDAL still reads what Janus does not cover (a file on the network, another kind of VRT), and the VRTs stay files GDAL and QGIS open.

11References

• Montero, D., Aybar, C., Mahecha, M. D., Martinuzzi, F., Söchting, M., & Wieneke, S. (2023). A standardized catalogue of spectral indices to advance the use of remote sensing in Earth system research. Scientific Data, 10, 197. The catalogue: awesome-spectral-indices (MIT); in Python, spyndex.
• Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83(1–2), 195–213.
• Hoyer, S., & Hamman, J. (2017). xarray: N-D labeled arrays and datasets in Python. Journal of Open Research Software, 5(1), 10.
• GDAL: VRT derived bands and pixel functions, gdal.org/drivers/raster/vrt.