plMapcalc — How to use

plMapcalc is a command-line tool. You specify input layers, an output layer, and a Python expression or script file. The expression is executed for every pixel in the raster, with access to the input values, output buffers, memory cells, and geographic coordinates.

Basic command syntax

shell shell

python plmapcalc.py [options] -i input1 [-i input2 …] [-o output1 …] -e 'expression'

The order of -i arguments determines the index used in the expression: the first -i is IN[0], the second is IN[1], etc.

Input layer specifier

syntax

-i filename.tif[:band[:shift_x[:shift_y]]]

Part	Default	Description
`filename.tif`	—	Path to any rasterio-readable file (GeoTIFF, VRT, …)
`band`	1	1-based band number to read
`shift_x`	0	Horizontal pixel offset. Negative = cells to the left of current cell.
`shift_y`	0	Vertical pixel offset. Negative = cells above current cell.

Layer shifting is used for neighbourhood operators (e.g. Laplace, slope). See Example 10.

Output layer specifier

syntax

-o filename.tif[:dtype[:nodata[:compress]]]

Part	Default	Values
`dtype`	float64	`uint8` `int16` `uint16` `int32` `uint32` `float32` `float64`
`nodata`	none	Any numeric value, e.g. `-9999` or `0`
`compress`	none	`DEFLATE` `DEFLATE2` `DEFLATE3` `LZW` `LZW2` `LZW3`

All command-line options

Short	Long	Description
`-i`	`--input`	Input layer (repeatable). See specifier above.
`-o`	`--output`	Output layer (repeatable). See specifier above.
`-e`	`--execute`	Python expression executed per pixel.
`-p`	`--program`	Path to a `.mc` or `.py` file with cell code.
	`--execute-begin`	Python code executed once before scanning starts — outside the RESTART loop. Use for MEM[] initialisation and one-time setup. Not called again on RESTART.
	`--execute-begin`	`--program-begin` — path to a `.mc` file with begin code.
	`--program-begin`	File with begin code.
	`--execute-end`	Python code executed once after each full scan. Call `RESTART()` here to trigger another pass. Called on every pass, unlike Begin.
	`--program-end`	File with end code.
`-m`	`--memory`	Number of memory cells to allocate (float64 array MEM).
`-r`	`--memory-read`	Load memory cells from a text file before processing.
`-s`	`--memory-store`	Save memory cells to a text file after processing.
`-0`	`--store-zeroes`	Include zero-valued cells in the `-s` output.
`-t`	`--threads`	Number of CPU threads for parallel JIT kernel (default: 1). Expressions that write to `MEM[]` always run single-threaded regardless of this setting — concurrent writes to shared memory cells are unsafe.
`-f`	`--force`	Overwrite existing output files.
`-q`	`--quiet`	Silent mode — suppresses all output including progress. Use twice (`-q -q`) for progress-only mode (no messages, progress bar visible).
	`--use-nan`	Treat NaN as a valid value instead of nodata.
	`--ignore-georef`	Skip projection and extent checks between input layers.
	`--crs`	Target CRS for all input layers, e.g. `EPSG:4326`. Layers with a different CRS are reprojected on-the-fly via WarpedVRT. Warning: reprojection is ~3–4× slower — pre-reproject with `gdalwarp` for large datasets.
	`--resolution`	Target pixel size: `ref` (default, first layer), `highest`, `lowest`, or a numeric value in CRS units. Activates WarpedVRT alignment. Warning: `highest` increases computation as the square of the scale factor.
	`--extent`	Output spatial extent: `ref` (default, first layer), `intersection`, or `union` (missing areas filled with nodata). Activates WarpedVRT alignment.
	`--resample`	Resampling method for WarpedVRT: `nearest` (default, recommended for discrete data), `bilinear`, `cubic`, `average`, `mode`, `lanczos`.
	`--gpu`	Run the pixel kernel on a CUDA GPU via `numba.cuda`. Falls back to CPU automatically if CUDA is unavailable. Expressions that write to `MEM[]` (e.g. `MEM[0] += IN[0]`) fall back to serial CPU — read-only access to `MEM[]` works on GPU.
	`--no-jit`	Disable numba JIT compilation and CUDA. The cell expression runs as plain Python — full Python syntax and all numpy functions are available, including dynamic slices (`IN[i:i+3]`), `dict`, `print()`, etc. Slower than JIT but useful for complex expressions or debugging. Implies `--threads 1`.

Variables available in expressions

Name	Type	Description
`IN[i]`	float64	Value of i-th input layer at current pixel. Read-only.
`OUT[i]`	float64	Output value for i-th output layer. Write to this.
`MEM[i]`	float64	Persistent memory cell (shared across all pixels and passes).
`GEOTRANS[i]`	float64	GDAL-style geotransform: [x_ul, cell_w, 0, y_ul, 0, cell_h].
`COL`	int	0-based column index of current pixel.
`ROW`	int	0-based row index of current pixel.
`COLS`	int	Total number of columns in the raster.
`ROWS`	int	Total number of rows in the raster.
`INPNUM`	int	Number of input layers.
`OUTNUM`	int	Number of output layers.
`MEMNUM`	int	Number of memory cells.

Built-in functions

Function	Available in	Description
`RESTART()`	execute-end	Trigger another full scan from row 0. Increments the iteration counter.
`ITERATION()`	all phases	Returns the current pass number (starts at 1).

Standard Python modules math and numpy (np) are pre-imported in all code phases.

Expression language

Expressions are valid Python. They run in a @njit-compiled context (when numba is available), which supports:

Arithmetic: + - * / // % **
Comparison and boolean: == != < > and or not
Math functions via math.sqrt, math.floor, math.log, etc.
if / elif / else
for loops with range()
Local variable assignment
Array indexing: IN[i], MEM[i]
Local Python lists — created and used inside the expression
G — shared global variable proxy (see below)

cell expression — local list python

values = []
for i in range(INPNUM):
    values.append(IN[i])
OUT[0] = sum(values) / INPNUM

Shared global variables — the `G` proxy

Variables assigned inside exec_cell are local to that call and are not visible in --execute-begin or --execute-end. Use the built-in G proxy object to share state across all three phases: G.name = value sets a shared variable, G.name reads it. Works identically with -e and -p.

--execute-begin

G.total = 0.0
G.count = 0

--execute (cell, called per pixel)

G.total += IN[0]
G.count += 1

--execute-end

MEM[0] = G.total / G.count if G.count > 0 else 0.0

ⓘ

MEM[] cells are also shared across all phases — they persist between passes (RESTART) and can be saved/loaded with -s/-r. MEM[] is NOT automatically zeroed between passes. Initialise MEM[] in --execute-begin — it runs only once before any scanning, so there is no need to guard with ITERATION()==1:

--execute-begin (called once, before pass 1)

MEM[0] = 0          # counter — stays intact between passes
MEM[1] =  1e308   # min — updated in cell, read in end
MEM[2] = -1e308   # max

Use G for lightweight intra-run state, MEM[] for persistent or multi-pass state.

⚠

G is not available in numba JIT mode. The G proxy uses a Python dict internally, which JIT cannot compile. ITERATION() and RESTART() are available in JIT mode. ITERATION() returns the correct pass number — it is passed as a parameter to the compiled kernel by the Python wrapper. RESTART() is a no-op in the JIT kernel itself — the restart signal is detected by the Python wrapper from the kernel return value. GPU mode (--gpu) does not support ITERATION()-based branching — the GPU stub always returns 1. Use CPU JIT (default) or --no-jit for multi-pass expressions. If your cell expression uses G, add --no-jit or replace G.var with MEM[i].

ℹ

IN[], OUT[], MEM[] are numpy arrays — they have no .append() method (this is a Python/numpy rule, not a numba restriction). Use index assignment: OUT[0] = value.
Local lists created inside the expression work fine in JIT — they are typed by numba automatically as long as they contain a single numeric type.

⚠

When using numba JIT, only numba-supported Python constructs work in the cell expression (i.e. the code passed to -e or -p). Not supported in JIT: dict, set, print(), string formatting, external function calls.
The --execute-begin and --execute-end phases run as regular Python and support the full language.
Use --no-jit to run the cell expression as plain Python when you need full language support — at the cost of slower processing.

⚠

Dynamic slicing inside loops — np.sum(IN[i:i+3]), np.mean(IN[i:i+3]) and similar calls with variable slice bounds inside a for loop are unreliable in numba JIT (behaviour varies between numba versions). Replace them with explicit element-wise arithmetic (see examples below). Also avoid np.array(existing_array) inside JIT — np.zeros(n) already returns an ndarray and needs no conversion. Alternatively, use --no-jit to run the expression as plain Python — dynamic slices and all numpy functions work without restriction, at the cost of slower processing.

instead of np.sum(IN[i:i+3])

val = IN[i] + IN[i+1] + IN[i+2]

instead of np.mean(IN[i:i+3])

val = (IN[i] + IN[i+1] + IN[i+2]) / 3.0

Automatic spatial alignment (`--resolution`, `--extent`, `--crs`)

By default plmapcalc requires all input layers to have identical CRS, pixel size, and extent. When any of --resolution, --extent, or --crs is specified, mismatched layers are aligned on-the-fly using rasterio.vrt.WarpedVRT. No temporary files are created — resampling happens lazily as each tile is read.

⚠

CRS reprojection is significantly slower. When --crs causes any layer to be reprojected, expect ~3–4× more I/O time per tile. For large datasets, pre-reproject your inputs with gdalwarp and run plmapcalc without --crs. plmapcalc will warn you when reprojection is active.

⚠

--resolution highest multiplies computation. If your inputs have resolutions of 10 m and 1 km, the output grid will be 100× wider and taller than the 1 km layer — 10 000× more pixels to compute. The default (--resolution ref) uses the first layer’s resolution.

ⓘ

For discrete data (class maps, patch IDs, hexagonal grids) always use --resample nearest (the default). Interpolating between integer IDs produces nonsensical intermediate values.

ⓘ

MEM[] auto-sizing: if -m is not set but the expression contains literal indices like MEM[4008], plmapcalc automatically allocates max_index + 1 cells and prints: [INFO] MEM[] auto-sized to 4009 cells. Use -m N explicitly when the index is dynamic (e.g. MEM[COLS * ROW + COL]).

Memory file format

Files used with -r (read) and -s (store) are plain text. Each line contains an index and a value separated by whitespace:

memory file

0  1.0
1  2.5
5  -9999.0
42 3.14159265358979

Indices do not need to be contiguous. Memory cells not listed in the file are initialised to 0. The array size must be declared with -m N where N is at least max_index + 1.

No-data handling

When a pixel in any input layer contains the no-data value (as declared in the GeoTIFF header), the cell expression is not executed for that pixel. The output cells are filled with the output no-data value instead.

Use --use-nan to disable this behaviour and treat NaN as a regular value.

Multi-pass processing

Calling RESTART() inside --execute-end causes the tool to repeat the entire scan from row 0. ITERATION() returns the current pass number so the cell expression can branch differently on each pass:

shell shell

python plmapcalc.py -i patches.tif -m 2000 -o areas.tif:float64:-9999:LZW \
  -e 'if ITERATION()==1: MEM[int(IN[0])] += 1.0
else: OUT[0] = MEM[int(IN[0])]' \
  --execute-end 'if ITERATION()==1: RESTART()'

Using macro files (-p)

For complex expressions, store the code in a file and pass it with -p. The file may have any extension (.mc or .py are conventional). The entire file contents are used as the cell body.

hexagon.mc .mc

# Size of hexagon side in pixels
size_x = 100
size2  = size_x // 2
size_y = int(math.sqrt(3.0) * size_x / 2.0)

col = COL + 1
row = ROW + 1
R = row // size_y
C = col // size_x
Cols = COLS // size_x + 1

R_center = R * size_y + size2
C_center = C * size_x
min_d = 2 * size_x * size_x
OUT[0] = 0.0

for r in range(-1, 2):
    Cc = C_center + ((R + r) % 2) * size2
    Rc = R_center + r * size_y
    if Cc < 0 or Cc >= COLS or Rc < 0 or Rc >= ROWS:
        continue
    dy = Rc - row
    for c in range(-1, 2):
        dx = Cc + c * size_x - col
        d  = dx*dx + dy*dy
        if d < min_d:
            min_d = d
            OUT[0] = 1 + C+c + (R+r)*Cols

shell shell

python plmapcalc.py -i reference.tif -o hexagons.tif:uint32:0:DEFLATE2 -p hexagon.mc

Windows notes

Use double quotes around -e expressions: -e "OUT[0] = IN[0] * 2"
Use backslash as path separator: -i C:\data\layer1.tif
The colon separator in layer specifiers still uses :: -i C:\data\layer1.tif:1:0:0

How to use

Basic command syntax

Input layer specifier

Output layer specifier

All command-line options

Variables available in expressions

Built-in functions

Expression language

Shared global variables — the G proxy

Automatic spatial alignment (--resolution, --extent, --crs)

Memory file format

No-data handling

Multi-pass processing

Using macro files (-p)

Windows notes

Shared global variables — the `G` proxy

Automatic spatial alignment (`--resolution`, `--extent`, `--crs`)