How do I analyse a Value Map in Python?

A Value Map is a Client File holding one attribute as a raster over a geographic area and a series of times, with a generation block on the record describing how to read it: what the numbers mean and in what units, the extent they cover, the tiles they are cut into, and where every raster sits in the blob. A simulation writes one per attribute, so a flood model might publish water level, velocity and wave height as separate Value Maps over the same area.

The question this answers is "what did the event touch". The Value Map says what the values were and where, and a Data Lab query says which of your Entities were there, so the two together say which records the event reached and when.

Reading one needs the optional extra, which a base install does not have.

pip install bruce-models[valuemap]

Finding a Value Map

Value Maps are Client Files, so they are found the same way as any other file. Each one names the attribute and the source it was generated from under Identity in its generation block, and those two together identify it.

import bruce_models as BM
# Not exported from the package root, so it is imported by path.
from bruce_models.client_file.client_file_value_map import ClientFileValueMap

# Either a session token from a login or a long-lived access token.
TOKEN = "YOUR_TOKEN"

bruce = BM.BruceApi({
    "account_id": "YOUR_ACCOUNT_ID",
    "session_id": TOKEN
})

# Every Value Map on the account, newest first.
query = (BM.ClientFileQueryBuilder(bruce)
         .purpose(BM.ClientFile.PURPOSE_VALUE_MAP_SEARCH)
         .newest_first())

for record in query.files():
    generation = (record.get("Data") or {}).get("generation") or {}
    identity = generation.get("Identity") or {}
    series = generation.get("Series") or {}
    print(f"{record.get('ID')}  "
          f"attribute={identity.get('Attribute')}  "
          f"dataset={identity.get('Dataset')}  "
          f"times={len(series.get('Times') or [])}")

# When you know the source and attribute, ask for that one directly.
one = BM.ClientFile.find_generated_value_map(
    bruce,
    entity_type_source_id=YOUR_SOURCE_ID,
    attribute="water_level"
)

Reading what it holds

Opening one downloads the blob once and reads the generation block off the record. From there everything is local, so walking every frame costs no requests at all.

Check is_complete before you read anything. It says whether every timestep that was asked for was actually written, and a series with a gap in the middle of the event otherwise looks exactly like a shorter series.

Cells that report the same value in every frame carry nothing a change can be read from, whether that is a permanent feature or an artefact of the model, and they are excluded for you.

value_map = ClientFileValueMap.open(bruce, VALUE_MAP_FILE_ID)

print(f"{len(value_map)} frames of {value_map.attribute} in {value_map.units}")
print(f"{value_map.grid.width} x {value_map.grid.height} texels")
print(f"values {value_map.value_range[0]:.3f} to {value_map.value_range[1]:.3f}")

if not value_map.is_complete:
    print(f"only {value_map.written_frames} of "
          f"{value_map.requested_frames} timesteps were written")

# Stated by you, since what counts as present depends on the attribute.
THRESHOLD = 0.0

# The frames worth reading: the event arriving, its worst, and it passing.
chosen = value_map.frames_of_interest(THRESHOLD)
print("frames of interest:", list(chosen))

# Which texels one frame holds a value over, and what was in them.
over = value_map.mask_above(THRESHOLD, frame=chosen[0])
values, covered = value_map.frame(chosen[0])
print(f"{int(over.sum()):,} texels over the threshold")

# Texels holding a value in every frame, which is what never receded.
print(f"{int(value_map.always_above(THRESHOLD).sum()):,} always over it")

Your own numbers, where they were corrected

Values are delivered in whatever reference the platform draws them against, which is not always the one your data was authored in. Where a correction was applied the Value Map carries it per texel, so the reading you supplied is recoverable exactly rather than described in prose.

value_source says whether anything was applied.frame(index, source=True) takes it back off, and the same flag runs through mask_above, peak_at andmean_at, so a threshold stated against your own data is compared against your own data. On a Value Map that declares no correction the two are the same numbers.

print(value_map.vertical_datum, value_map.transform, value_map.source_range)

delivered, covered = value_map.frame(0)
mine, mine_covered = value_map.frame(0, source=True)

# A level you know from your own data, measured against your own data.
reached = value_map.mask_above(1.25, frame=0, source=True)

Sampling one place

When the question is about a site rather than a whole frame, sample it directly.sample_at reads one position andseries_at reads that position through the whole series, each returning the delivered number and your own alongside each other, so a script can report either or both without deciding up front which one it works in.

Both read only the tile the position falls in, so watching one site through a long series costs a fraction of walking the frames. sample_of andseries_of do the same for an Entity, which covers several texels and so reports a peak and a mean rather than a single reading.

A position outside the raster comes back as nothing at all, while one inside it that was never sampled comes back with covered false and no value. Those are different answers and it is worth reporting them differently.

reading = value_map.sample_at(174.76, -36.85, frame=0)
if reading is None:
    print("outside the modelled area")
elif not reading["covered"]:
    print("inside it, but nothing was modelled there")
else:
    print(f"{reading['time']}: {reading['value']:.3f} {reading['units']} "
          f"(yours: {reading['source_value']:.3f})")

# The same place through every timestep.
series = value_map.series_at(174.76, -36.85)
worst = max((r for r in series if r["value"] is not None),
            key=lambda r: r["value"])
print(f"peaked at {worst['value']:.3f} on {worst['time']}")

# One record, which covers several texels.
summary = value_map.sample_of({"Bruce": item})
print(f"{summary['covered']} of {summary['texels']} texels modelled, "
      f"peak {summary['peak']}, mean {summary['mean']}, "
      f"peak in your own numbers {summary['source_peak']}")

What the event reached

An Entity resolves to the texels it covers, and the raster does not move, so list the Entities once, index them by their texels, and every frame after that is answered locally. The cost then follows the size of the event rather than the number of records, and the match is exact against the raster.

shapes_above is the other way round: it traces each area of a frame into a Bruce geometry you can send to Data Lab as a spatial criterion. Reach for it when you want the outline itself, to draw or to hand to something else.

import numpy

PARCEL_TYPE_ID = "YOUR_ENTITY_TYPE_ID"

parcels = BM.DataLabBuilder(bruce, page_size=1000).entity_type(PARCEL_TYPE_ID)
timestamps = value_map.timestamps()

# Listed once. The Bruce select is enough, since only the identity and
# where it sits are needed.
records, by_texel = [], {}
for page in parcels.pages("Bruce"):
    for item in page.items:
        texels = value_map.texels_of({"Bruce": item})
        if not texels:
            continue
        position = len(records)
        records.append(({"Bruce": item}, texels))
        for texel in texels:
            by_texel.setdefault(texel, []).append(position)

print(f"{len(records):,} records on the raster")

# Each reading is its own edit, so a record keeps one per frame.
changes = BM.ChangeSetBuilder(bruce, allow_multi_same_entity=True)

for index in chosen:
    when = timestamps[index]
    values, covered = value_map.frame(index)
    over = value_map.mask_above(THRESHOLD, frame=index)

    # Walked from the texels over the threshold rather than over every record.
    reached = set()
    rows, columns = numpy.nonzero(over)
    for row, column in zip(rows.tolist(), columns.tolist()):
        reached.update(by_texel.get((row, column), ()))

    for position in reached:
        entity, texels = records[position]
        level = value_map.peak_at(texels, values=values, covered=covered)
        if level is None:
            continue
        reading = {}
        BM.Entity.set_value(reading, "flood/time", when)
        BM.Entity.set_value(reading, "flood/level", level)
        changes.queue(entity, reading)

    print(f"frame {index} ({when}): {len(changes):,} readings")

print(f"{len(changes.entity_ids):,} Entities reached")

# Written as ordinary updates, since a saved Change Set holds one item
# per Entity and these are a series per record.
report = changes.apply()
print(f"{report['succeeded']:,} readings saved")

Things worth knowing

Nothing states what the numbers mean, so the threshold is yours to name. The Value Map publishes what its numbers are worth and in what units, not whether a value counts as present. Zero is the line for something measured up from whatever is beneath it, such as a depth, and is merely wherever the datum sits for a height. Leave the threshold out and one is read off the distribution instead, which separates whichever region dominates the range rather than the event you are looking for.

A texel that was never sampled is not a reading of zero. Every raster carries its own coverage, which is the second value frame returns, and a record over unsampled ground has nothing to report. Coverage belongs to the raster you read it from, so never mask one raster's numbers with another's coverage.

The outline can never be finer than a texel. Areas are traced off the raster, so the shape you query with steps around whole cells. A texel can be hundreds of metres across, which is much larger than a building or a parcel.

Values are quantised over the archive's own range. Depending on how wide that range is, a value is packed into either one channel or two, so the step is the range divided by 255 or by 65,535. Read the numbers through the library rather than off the pixels: the same picture read as one channel where it holds two is wrong by up to the whole range.

Water that never recedes is still water. A record standing in a permanent feature reads as covered in every frame, so a count below the length of the series is a record the event arrived at. The reverse is weaker: a full count only says some part of it was covered throughout, which is true both of a record wholly inside the feature and of one straddling its edge. Decide which of the two you are reporting before writing anything, because both are legitimate answers to different questions.

Repeated readings on one Entity need allow_multi_same_entity. A Change Set built without it merges edits per Entity, so a series of timestamped readings would collapse to the last one. A builder created with it can only be applied, not saved as a Change Set record, because a record holds one item per Entity.