plMapcalc — Examples

The following examples show different approaches to using plMapcalc, from the simplest single-expression calls to multi-pass algorithms with memory buffers. They are ordered from simplest to most complex.

All examples are written for Linux/macOS (single quotes around -e arguments). On Windows, use double quotes instead.

Contents

Adding two or more layers
Calculation of NDVI
Creating a mask file
Statistics calculation
Reclassification of input layer
Area calculation and assigning
Space partitioning
Thiessen (Dirichlet/Voronoi) tessellation
Image histogram matching and equalization
Laplace operator for edge detection
Boundary detection for discrete data (patch IDs, hexagonal grids)
Köppen-Geiger climate classification

Example 1Adding two or more layers

Assume you have two layers: layer1.tif and layer2.tif with the same resolution, extent, and projection. In the case of integer data, plMapcalc reads them as float64 internally.

shell shell

python plmapcalc.py -i layer1.tif -i layer2.tif \
                  -o sum.tif \
                  -e 'OUT[0] = IN[0] + IN[1]'

The output layer is created with the default options: float64 data type, no no-data value. IN[0] refers to the first input layer, IN[1] to the second.

To write integer output with a no-data value:

shell shell

python plmapcalc.py -i layer1.tif -i layer2.tif \
                  -o sum.tif:int32:-9999 \
                  -e 'OUT[0] = IN[0] + IN[1]'

To control rounding explicitly (default truncates toward zero):

shell shell

python plmapcalc.py -i layer1.tif -i layer2.tif \
                  -o sum.tif:int32:-9999 \
                  -e 'OUT[0] = math.floor(IN[0] + IN[1])'

To sum an arbitrary number of layers use INPNUM:

shell shell

python plmapcalc.py -i layer1.tif -i layer2.tif -i layer3.tif \
                  -o sum.tif:float32:-9999 \
                  -e '
for i in range(INPNUM):
    OUT[0] += IN[i]
'

Example 2Calculation of NDVI

The Normalised Difference Vegetation Index from Landsat bands 4 and 3:

NDVI = (band4 − band3) / (band4 + band3)

The denominator can be zero, which would produce ±Inf. Handle the edge case explicitly:

shell shell

python plmapcalc.py -i band4.tif \
                  -i band3.tif \
                  -o ndvi.tif:float32:-9999 \
                  -e '
x = IN[0] + IN[1]
OUT[0] = -9999 if x == 0.0 else (IN[0] - IN[1]) / x
'

The local variable x is private to the current pixel. If x is zero the output is set to the no-data value -9999.

Example 3Creating a mask file

Create a binary mask where pixels equal to a specific class value (e.g. country index 38) are 1 and all others are 0:

shell shell

python plmapcalc.py -i europe.tif \
                  -o my_country_mask.tif:uint8 \
                  -e 'OUT[0] = 1 if IN[0] == 38 else 0'

The result contains 1 over the target country and 0 elsewhere.

Example 4Statistics calculation

Use the NDVI layer from Example 2 and the mask from Example 3 to compute the average NDVI for the selected country. The memory buffer accumulates the count and sum while scanning, and the average is computed in --execute-end:

shell shell

python plmapcalc.py -i ndvi.tif -i my_country_mask.tif \
                  -m 3 \
                  -e '
if IN[1] > 0:
    MEM[0] += 1.0
    MEM[1] += IN[0]
' \
                  --execute-end 'MEM[2] = MEM[1] / MEM[0] if MEM[0] > 0 else 0' \
                  -s result.txt

The saved file result.txt contains three lines: MEM[0] — pixel count, MEM[1] — sum of values, MEM[2] — average. There is no raster output in this call.

💡

When a pixel in any input layer has the no-data value, the cell expression is skipped entirely. So ndvi.tif no-data pixels are already excluded without an explicit check.

Example 5Reclassification of input layer

Reclassify a 16-class land-cover map to 6 general classes using a lookup table stored in a memory file. The input classes are:

11 – Open Water, 12 – Perennial Ice/Snow
21–24 – Developed (various intensities)
31 – Barren Land, 41–43 – Forest types
52 – Shrub, 71 – Grassland, 81–82 – Agriculture
90, 95 – Wetlands

Create mem_map.txt with the mapping (index → output class):

mem_map.txt memory file

shell shell

python plmapcalc.py -i lc_classes.tif \
                  -m 100 -r mem_map.txt \
                  -o lc_reclass.tif:uint8:0 \
                  -e 'OUT[0] = MEM[int(IN[0])]'

The expression int(IN[0]) converts the float pixel value to an integer index into MEM[]. MEM cells not in the file default to 0 (no-data).

Figure 1. Reclassification result: original 16-class land cover (left) and reclassified 6-class map (right).

Example 6Area calculation and assigning

Count how many pixels belong to each patch, then assign the area value back to every pixel in each patch. Two-step approach:

Step 1 — count pixels per patch

shell shell

python plmapcalc.py -i patches.tif \
                  -m 2000 \
                  -e 'MEM[int(IN[0])] += 1.0' \
                  -s areas.txt

Step 2 — assign areas

shell shell

python plmapcalc.py -i patches.tif \
                  -r areas.txt -m 2000 \
                  -o areas.tif:float64:-9999:LZW \
                  -e 'OUT[0] = MEM[int(IN[0])]'

Both steps can be combined in a single call using ITERATION() and RESTART():

shell — one-call variant 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()
'

Figure 2. Area calculation result for forest districts in Poland. Colour encodes patch area: white = smallest, dark red = largest.

Example 7Space partitioning

Assign a unique patch ID to each pixel based on its position. The rook-topology square grid divides the raster into tiles of size × size pixels:

shell shell

python plmapcalc.py -i reference.tif \
                  -o sq_patches.tif:uint32:0:DEFLATE \
                  -e '
size = 100
OUT[0] = 1 + COL // size + (ROW // size) * (1 + COLS // size)
'

For a hexagonal grid, store the more complex expression in a file:

hexagon.mc .mc

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

Figure 3. Square patches with rook topology. Patches are filled with random colours; labels show patch IDs.

Figure 4. Hexagonal patches generated by hexagon.mc. Patches are filled with random colours; labels show patch IDs.

Example 8Thiessen (Dirichlet/Voronoi) tessellation

Create a Voronoi diagram from a set of seed points. First, prepare a memory file with alternating X, Y coordinates (one pair per seed). Save it as xy_coordinates.txt:

xy_coordinates.txt memory file

 0 4749361.26105342
 1 3114579.65054740
 2 4759413.25837782
 3 3123209.19542023
 4 4764059.93638626
 5 3098648.18308986
 6 4751068.20399530
 7 3097036.07031142
 8 4778568.95139223
 9 3111924.40597114
10 4789948.57100476
11 3096182.59884048
12 4781129.36580505
13 3066121.43703073
14 4793077.96639820
15 3073518.18977887
16 4776862.00845036
17 3088596.18576547
18 4790991.70280257
19 3121976.40329554

Even indices hold X (easting), odd indices hold Y (northing). The --execute-begin phase converts geographic coordinates to pixel column/row; the cell expression then finds the nearest seed:

shell shell

python plmapcalc.py -i reference.tif \
                  -o thiessen.tif:uint32:0:DEFLATE2 \
                  -r xy_coordinates.txt \
                  --execute-begin '
i = 0
while i < MEMNUM:
    x = MEM[i]
    y = MEM[i+1]
    MEM[i]   = (x - GEOTRANS[0]) / GEOTRANS[1]
    MEM[i+1] = (y - GEOTRANS[3]) / GEOTRANS[5]
    i += 2
' \
                  -e '
min_d = 2.0 * ROWS * COLS
OUT[0] = 0.0
i = 0
while i < MEMNUM:
    dy = MEM[i+1] - ROW
    dx = MEM[i]   - COL
    d  = dx*dx + dy*dy
    if d < min_d:
        min_d = d
        OUT[0] = 1 + i // 2
    i += 2
'

ℹ

This example writes to MEM in --execute-begin, not in the cell expression, so the parallel JIT kernel is still used for the pixel loop.

Figure 5. Thiessen (Voronoi) tessellation. Polygons are filled with random colours; labels show polygon IDs. Seed points are marked in grey.

Example 9Image histogram matching and equalization

Histogram equalization stretches the contrast of a greyscale image by mapping pixel values through the image's cumulative distribution function (CDF). A two-pass approach is used: pass 1 builds the histogram; pass 2 applies the reclassification.

shell — histogram equalization shell

python plmapcalc.py --input=greyscale_input.tif \
                  --memory=1000 \
                  --output=grey_eq_output.tif:int32 \
                  --execute-begin '
for i in range(256):
    MEM[i] = i / 256.0
' \
                  --execute '
if ITERATION() == 1:
    MEM[int(IN[0]) + 500] += 1.0
else:
    OUT[0] = MEM[int(IN[0]) + 500]
' \
                  --execute-end '
if ITERATION() == 1:
    for i in range(1, 256):
        MEM[i + 500] += MEM[i - 1 + 500]
    v = MEM[255 + 500]
    for i in range(256):
        MEM[i + 500] /= v
    for i in range(500, 500 + 256):
        best_j, best_v = 0, abs(MEM[0] - MEM[i])
        for j in range(1, 256):
            a = abs(MEM[j] - MEM[i])
            if a < best_v:
                best_j, best_v = j, a
        MEM[i] = best_j
    RESTART()
'

For histogram matching (match one image's histogram to a reference image), add the reference as a second input and accumulate both histograms in pass 1:

shell — histogram matching shell

python plmapcalc.py --input=greyscale_input.tif \
                  --input=greyscale_reference.tif \
                  --memory=1000 \
                  --output=grey_matched_output.tif:int32 \
                  --execute-begin '
for i in range(256):
    MEM[i] = i / 256.0
' \
                  --execute '
if ITERATION() == 1:
    MEM[int(IN[0]) + 500] += 1.0
    MEM[int(IN[1])]       += 1.0
else:
    OUT[0] = MEM[int(IN[0]) + 500]
' \
                  --execute-end '
if ITERATION() == 1:
    for i in range(1, 256):
        MEM[i + 500] += MEM[i - 1 + 500]
        MEM[i]       += MEM[i - 1]
    v_in  = MEM[255 + 500]
    v_ref = MEM[255]
    for i in range(256):
        MEM[i + 500] /= v_in
        MEM[i]       /= v_ref
    for i in range(500, 500 + 256):
        best_j, best_v = 0, abs(MEM[0] - MEM[i])
        for j in range(1, 256):
            a = abs(MEM[j] - MEM[i])
            if a < best_v:
                best_j, best_v = j, a
        MEM[i] = best_j
    RESTART()
'

Figure 6. Histogram equalization: original greyscale image (left) and equalized output (right) with stretched contrast.

Figure 7. Histogram matching result. Left: reference image. Centre: input image before matching. Right: input image after its histogram has been matched to the reference.

Example 10Laplace operator for edge detection

The discrete 8-direction Laplace kernel:

kernel

-1  -1  -1
-1   8  -1
-1  -1  -1

plMapcalc implements neighbourhood operations through layer shifting. The same file is passed nine times with different shift_x:shift_y offsets. A negative shift means "take from cells to the left / above":

shell shell

python plmapcalc.py \
  -i input.tif \
  -i input.tif::-1:-1 -i input.tif::-1:0 -i input.tif::-1:1 \
  -i input.tif::0:-1                      -i input.tif::0:1 \
  -i input.tif::1:-1  -i input.tif::1:0  -i input.tif::1:1 \
  -o output.tif:uint8 \
  -e '
OUT[0] = 8 * IN[0]
for i in range(1, 9):
    OUT[0] -= IN[i]
'

IN[0] is the centre pixel. IN[1]–IN[8] are the eight neighbours at the offsets listed. The band index is omitted (defaults to band 1).

ℹ

The Laplace kernel is designed for continuous data (elevation, imagery). It detects smooth gradients, not sharp ID boundaries. For discrete data (patch IDs, class maps, hexagonal grids) use the boundary detection method below.

Edge-detected result after Laplace operator

Figure 8. Laplace operator for edge detection. Left: source greyscale image. Right: edges extracted by the 8-direction Laplace kernel.

Example 11Boundary detection for discrete data (patch IDs, hexagonal grids)

The Laplace kernel does not produce correct boundaries for discrete ID rasters (patch maps, land cover classes, hexagonal grids) because it detects gradients, not value changes. A pixel lies on a boundary when its value differs from any of its neighbours:

shell shell

python plmapcalc.py \
  -i patches.tif \
  -i patches.tif::-1:-1 -i patches.tif::-1:0 -i patches.tif::-1:1 \
  -i patches.tif::0:-1                         -i patches.tif::0:1 \
  -i patches.tif::1:-1  -i patches.tif::1:0   -i patches.tif::1:1 \
  -o edges.tif:uint8:255:DEFLATE \
  -e '
is_edge = 0
for i in range(1, 9):
    if IN[i] != IN[0]:
        is_edge = 1
        break
OUT[0] = is_edge
'

Output pixels equal 1 on patch boundaries and 0 inside patches. This correctly detects all edges — including vertical boundaries — on hexagonal grids, square patch maps, and any discrete classification raster.

💡

For hexagonal grids the centre pixel physically belongs to one hexagon while its diagonal neighbours may belong to different hexagons. The 8-neighbour comparison above covers all six hexagon sides because the hex grid geometry guarantees that at least one of the 8 raster neighbours crosses each side boundary.

Example 11Köppen-Geiger climate classification

The Köppen-Geiger classification requires monthly temperature and precipitation for all 12 months. Download the data from WorldClim and store as:

Temperature: t01.tif … t12.tif (in tenths of °C)
Precipitation: p01.tif … p12.tif (mm/month)

The classification logic (Spinoni et al. 2015) is too long for an inline -e argument, so it is stored in a macro file. Save the following as macros/KG_classification_13.mc:

KG_classification_13.mc .mc

# KG class constants
EF,ET,BW,BS,Am,Aw,Af = 1,2,3,4,5,6,7
CS,CW,CF,DS,DW,DF    = 8,9,10,11,12,13

Tcold = 10000.0; Thot = -10000.0; Tmon10 = 0.0
MAT = 0.0; MATw = 0.0; MATs = 0.0

for i in range(12):
    Ta = IN[i]
    if Ta < Tcold: Tcold = Ta
    if Ta > Thot:  Thot  = Ta
    if Ta > 100.0: Tmon10 += 1.0
    MAT += Ta
    if 3 < i+1 < 10: MATs += Ta
    else:           MATw += Ta

Tcold  *= 0.1; Thot *= 0.1; MAT /= 120.0

MAP = 0.0; MAPw = 0.0; MAPs = 0.0
Pdry = 10000.0; Pwdry = 10000.0; Psdry = 10000.0
Pwwet = 0.0;  Pswet = 0.0

for i in range(12):
    Pr = IN[i + 12]
    if Pr < Pdry: Pdry = Pr
    MAP += Pr
    if 3 < i+1 < 10:
        if Pr < Psdry: Psdry = Pr
        if Pr > Pswet: Pswet = Pr
        MAPs += Pr
    else:
        if Pr < Pwdry: Pwdry = Pr
        if Pr > Pwwet: Pwwet = Pr
        MAPw += Pr

if MATw > MATs:
    MAPw,MAPs   = MAPs,MAPw
    Pwdry,Psdry = Psdry,Pwdry
    Pwwet,Pswet = Pswet,Pwwet

if   MAPw > 0.7*MAP: Pthre = 20.0*MAT
elif MAPs > 0.7*MAP: Pthre = 20.0*MAT + 280.0
else:               Pthre = 20.0*MAT + 140.0

alpha = (Psdry < 30.0) and (Psdry < Pwwet/3.0) and (MAPs < MAPw)
beta  = (Pwdry < Pswet/10.0) and (MAPw < MAPs)
delta = not (alpha and beta)

KG = 0
if Thot < 10.0:
    KG = ET if Thot > 0.0 else EF
elif Pthre >= MAP:
    KG = BW if MAP < Pthre/2.0 else BS
elif Tcold >= 18.0:
    if Pdry < 60.0:
        KG = Am if Pdry >= 100.0 - MAP/25.0 else Aw
    else: KG = Af
elif Tcold < -3.0:
    KG = DS if alpha else (DW if beta else (DF if delta else 0))
else:
    KG = CS if alpha else (CW if beta else (CF if delta else 0))

OUT[0] = float(KG)

shell shell

python plmapcalc.py \
  -i t01.tif -i t02.tif -i t03.tif -i t04.tif -i t05.tif -i t06.tif \
  -i t07.tif -i t08.tif -i t09.tif -i t10.tif -i t11.tif -i t12.tif \
  -i p01.tif -i p02.tif -i p03.tif -i p04.tif -i p05.tif -i p06.tif \
  -i p07.tif -i p08.tif -i p09.tif -i p10.tif -i p11.tif -i p12.tif \
  -o kg13.tif:uint8:0:LZW \
  -p macros/KG_classification_13.mc

24 input files are read simultaneously; one output raster is created. The output is uint8 with 0 as no-data, compressed with LZW. The macro file is read from the macros/ subdirectory. The order of input files is significant: months 1–12 are temperature (IN[0]–IN[11]), months 1–12 are precipitation (IN[12]–IN[23]).

Köppen-Geiger climate classification map of the world produced by plMapcalc

Figure 9. Köppen-Geiger climate classification (13 classes) produced by plMapcalc from WorldClim monthly temperature and precipitation data. (Netzel & Stepinski, Journal of Climate, 2016.)