arrow/parquet as a metadata format

[sjperkins]

Following on from @bmerry's #4 (comment), I thought I'd spitball something arrow/parquet based for a couple of reasons:

1. Internally, parquet chunks data by row groups.
2. This allows readers to perform partial column reads

The following script spitballs placing MeerKAT sensor data into arrow tables an parquet files. It assumes that each row is associated with a dump -- perhaps this needs to be corrected? However, proceeding from this assumption, one can store baseline, channel or any other dimensionality data in each row.

What happens if data doesn't vary by row?: Dump it into the arrow table metadata as json/msgpack.

from itertools import product
from datetime import datetime
import json
from tempfile import TemporaryDirectory

import numpy as np
import pyarrow as pa
import pyarrow.compute as pac
import pyarrow.parquet as pq

ANTS = [f"m{x:03}" for x in range(7)]
POLS = ["hh", "hv", "vh", "vv"]
ANT1, ANT2 = np.triu_indices(len(ANTS), 0)
CORRPRODS = [(f"{ANTS[a1]}{p[0]}", f"{ANTS[a2]}{p[1]}") for (a1, a2), p in product(zip(ANT1, ANT2), POLS)]
START_TIME = datetime.now().timestamp()
NBL = len(ANT1)
NDUMPS = 100
DUMP_RATE = 8.0
NCHAN = 16

if __name__ == "__main__":
  time = START_TIME + np.arange(NDUMPS) * DUMP_RATE
  chan_freq = np.linspace(.856e9, 2 * .856e9, NCHAN)

  # timestamps are the table's primary row dimension
  telstate = pa.Table.from_pydict({
    "time": time,
    # Categorical data
    "state": pa.DictionaryArray.from_arrays(indices=np.random.randint(0, 4, NDUMPS), dictionary=pa.array(["slew", "track", "scan", "think"])),
    # Parquet can't write categorical of FSLA yet
    # "chan_freq": pa.DictionaryArray.from_arrays(indices=[0] * NDUMPS, dictionary=pa.FixedSizeListArray.from_arrays(chan_freq, NCHAN)),
    # Arrow can't create categorical of FSTA yet
    # "chan_freq": pa.DictionaryArray.from_arrays(indices=[0] * NDUMPS, dictionary=pa.FixedShapeTensorArray.from_numpy_ndarray(chan_freq.reshape((1, NCHAN)))),
    "m003_humidity": pa.FixedShapeTensorArray.from_numpy_ndarray(np.random.random((NDUMPS, NBL, NCHAN))),
    "m003_pointing_azel": pa.FixedShapeTensorArray.from_numpy_ndarray(np.random.random((NDUMPS, len(ANTS), NCHAN, 2)))
  },
  # Everything without a primary time dimension goes into the metadata section
  # This should probably be msgpack'd to avoid losing precision
  # If data is split across multiple parquet files, the metadata is replicated
  # in each, which may be undesirable
  metadata={
    "corrprods": json.dumps(CORRPRODS),
    "dump_rate": json.dumps(DUMP_RATE),
    "antenna": json.dumps(ANTS),
    "chan_freq": json.dumps(chan_freq.tolist()),
    "polarizations": json.dumps(POLS)
  })

  with TemporaryDirectory(prefix="telstate") as tmpdir:
    # Write out a single parquet file, in chunks of 10 rows
    pq.write_to_dataset(telstate, tmpdir, row_group_size=10)

    # Lazy view over the data
    ds = pq.ParquetDataset(
      tmpdir,
      filters=(
        (pac.field("time") > pac.scalar(time[3])) &
        (pac.field("time") < pac.scalar(time[-3]))
      )
    )

    # Perform projection, only read this column
    ds_table = ds.read(["m003_pointing_azel"])
    assert len(ds_table) == NDUMPS - 3 - 4

@ludwigschwardt Is there anything that wouldn't work in the above conceptualisation? I'd need to figure out if it's possible to remotely access a parquet file via http(s) (it's certainly possible to do so via S3 which is built on top of http(s) anyway.

[sjperkins]

I'd need to figure out if it's possible to remotely access a parquet file via http(s)

If it's not possible to do it with pyarrow's native handling, it has fsspec extension points, which allow the use of fsspec's http implementation.

[ludwigschwardt]

I don't know Parquet or Arrow at all but here are some thoughts:

• The native metadata format is telstate. The sensor timestamps are not aligned with the dumps or each other and part of the magic / burden of katdal is to make that happen...
• As much as possible is computed and not stored on purpose (cf. UVW and even timestamps...).
• If there will be one Parquet file per sensor, I would call it "unappetising" for a lack of a better word 😅 On the other hand, a single Parquet file would be grand!
• Some sensors have serious values and not just "floats" (e.g. arrays of complex numbers of shape (64, 2, 32768) in the case of BP cal solutions).
• One option could be to let katdal produce a "post-RDB" metadata representation with sensors aligned with dumps etc, a bit like the ASDM fillers of yore.

[ludwigschwardt]

Also, I'm getting the feeling that lazy loading of metadata is not trivial, especially if your file is on the network and you don't want to copy all of it across. Will the fsspec route only work if you have many small files? I see it has block-based caching for partial downloads too. If not trivial, I actually like Bruce's suggestion of a telstate-like interface over the network, resulting in a "lazy memory" RDB backend instead.

I think a guiding principle for me would be to see how an end user would actually use this metadata previously served by katdal and work back towards the data source before I can make up my mind about any possible intermediate formats.

[sjperkins]

• The native metadata format is telstate. The sensor timestamps are not aligned with the dumps or each other and part of the magic / burden of katdal is to make that happen...

Yep, I suspected this might be the case and if there isn't a "good" definition of what each row is, then it's probably not possible to store multiple sensors in one parquet file.

• If there will be one Parquet file per sensor, I would call it "unappetising" for a lack of a better word 😅 On the other hand, a single Parquet file would be grand!

I suspect to settle on a single Parquet file, one would have to have a uniform definition for what constitutes a row, and if sensors timestamps aren't aligned then multiple files per sensor would be needed, and I wouldn't find that appealing either.

• Some sensors have serious values and not just "floats" (e.g. arrays of complex numbers of shape (64, 2, 32768) in the case of BP cal solutions).

I do think this can be handled (see the pa.FixedShapeTensorArray above) even if the ndarray had an extra (2,) on the shape to represent the complex pair. It's all buffers under the hood anyway.

• One option could be to let katdal produce a "post-RDB" metadata representation with sensors aligned with dumps etc, a bit like the ASDM fillers of yore.

Could this happen in the file writer (if there was a parquet file writer)? Although this would expand the size of the stored metadata.

[sjperkins]

Also, I'm getting the feeling that lazy loading of metadata is not trivial, especially if your file is on the network and you don't want to copy all of it across. Will the fsspec route only work if you have many small files? I see it has block-based caching for partial downloads too. If not trivial, I actually like Bruce's suggestion of a telstate-like interface over the network, resulting in a "lazy memory" RDB backend instead.

All of the remote HTTP-based backends (HTTP/S3/GCS) are based on HTTP. AFAICT they read the parquet header to figure out the byte ranges of the columns and row groups, which are then passed in as a header to the HTTP GET requests. So it doesn't read the entire parquet file over HTTP

I think a guiding principle for me would be to see how an end user would actually use this metadata previously served by katdal and work back towards the data source before I can make up my mind about any possible intermediate formats.

I hadn't thought of presenting telstate to the end user specifically. I guess xarray-kat needs to consume telstate in a similar manner to katdal?

[bmerry]

I hadn't realised parquet didn't allow for multiple datasets in a single file. Perhaps something like sqlite or HDF5 is better suited then - although I don't know how well either of those will play with partial reads over HTTP.

[sjperkins]

I hadn't realised parquet didn't allow for multiple datasets in a single file.

Gemini suggested having (timestamp, sensor_id and metric_name) as the primary key:

Column Name	Data Type	Repetition Type	Arrow/Parquet Interpretation	Purpose
timestamp	TIMESTAMP	REQUIRED	Timestamp(unit='ns')	Exact measurement time.
sensor_id	STRING	REQUIRED	String	Unique sensor identifier.
metric_name	STRING	REQUIRED	String	What was measured (e.g., acceleration, temperature).
value_scalar	DOUBLE	OPTIONAL	Float64	For simple single-float readings (like temperature).
value_vec2d	LIST of DOUBLE	OPTIONAL	List<item: float64>	For 2D readings (e.g., Lat/Lon, or X/Y).
value_vec3d	LIST of DOUBLE	OPTIONAL	List<item: float64>	For 3D vectors (e.g., X/Y/Z acceleration).
value_vec4d	LIST of LIST of DOUBLE	OPTIONAL	List<item: List<item: float64>>	For 4d vectors..

Then, one of the value* columns would contain the actual data. As these columns can have null entries, this is apparently quite compact. @ludwigschwardt What would the full range of value types be?

[bmerry]

That seems like a rather messy solution. It doesn't give a natural way to get all the data for one sensor - you have to select it out of a single table.

Telstate encodes values with msgpack with some extensions to allow for things like numpy arrays and distinguishing tuples from lists. I'd suggest that for maximum compatibility, any alternative file format should continue to store sensor values as msgpack blobs rather than trying to shoehorn things into an Arrow/HDF5/JSON/etc data model.

[sjperkins]

That seems like a rather messy solution.

Yeah, fair enough. However, if the actual data payload is msgpacked/encoded anyway, how about the following:

Column Name	Data Type	Repetition Type	Arrow/Parquet Interpretation	Purpose
timestamp	TIMESTAMP	REQUIRED	Timestamp(unit='ns')	Exact measurement time.
sensor_id	STRING	REQUIRED	String	Unique sensor identifier.
sensor_type	INT	REQUIRED	Integer	Type Identifier describing the payload
payload	BYTE_ARRAY	REQUIRED	Binary	Raw binary data

It doesn't give a natural way to get all the data for one sensor - you have to select it out of a single table.

Are you referring to the snippet below to retrieve a single sensor (although multiple sensors could be read via the ds.read command)?

    # Lazy view over the data
    ds = pq.ParquetDataset(
      tmpdir,
      filters=(
        (pac.field("time") > pac.scalar(time[3])) &
        (pac.field("time") < pac.scalar(time[-3]))
      )
    )

    # Perform projection, only read this column
    ds_table = ds.read(["m003_pointing_azel"])

While the above probably isn't al that nice for a telstate consumer used to pulling data out of a Mapping interface, I do feel that something like the above could be encapsulated in a __getitem__ call. It should also only query the rows specified by the filters expression and the columns specified in ds.read

[sjperkins]

It should also only query the rows specified by the filters expression and the columns specified in ds.read

i.e. it should only pull this data out of a local parquet file, or over the network.

[sjperkins]

If there will be one Parquet file per sensor, I would call it "unappetising" for a lack of a better word 😅 On the other hand, a single Parquet file would be grand!

Let me play devils advocate here and ask why a parquet file per sensor would be a bad thing?

[sjperkins]

However, if the actual data payload is msgpacked/encoded anyway, how about the following:

Here's a more concrete example:

from datetime import datetime
from itertools import product
from random import randbytes, randrange
from uuid import uuid4
from tempfile import TemporaryDirectory

import msgpack
import numpy as np
import pyarrow as pa
import pyarrow.compute as pac
import pyarrow.parquet as pq

NSENSOR_SAMPLES = 100

if __name__ == "__main__":
  sensor_id = np.random.choice(["m003_humidity", "m003_pointing_azel"], NSENSOR_SAMPLES)
  timestamp = np.sort(np.random.random(size=NSENSOR_SAMPLES) * NSENSOR_SAMPLES)
  # Random payload bytes, but this would correspond to msgpack'd data
  payload = [msgpack.packb(uuid4().hex[:randrange(4, 10)]) for _ in range(NSENSOR_SAMPLES)]
  telstate = pa.Table.from_pydict({
    "sensor_id": sensor_id,
    "timestamp": timestamp,
    "payload": payload,
  })

  print(telstate)

  with TemporaryDirectory(prefix="telstate") as tmpdir:
    # Write out a single parquet file, in chunks of 10 rows
    pq.write_to_dataset(telstate, tmpdir, row_group_size=10)

    # Pull out data for the two sensors between these timestamps,
    ds = pq.ParquetDataset(
      tmpdir,
      filters=(
        (pac.field("sensor_id").isin(["m003_pointing_azel", "m003_humidity"])) &
        (pac.field("timestamp") > pac.scalar(timestamp[3])) &
        (pac.field("timestamp") < pac.scalar(timestamp[-3]))
      )
    )

    # Perform projection, only read this column
    ds_table = ds.read(["payload"])

[bmerry]

However, if the actual data payload is msgpacked/encoded anyway, how about the following:

Here's a more concrete example:

Isn't that going to have to read over the entire file in order to apply the filter, even if you're only interested in a subset of sensors (or a single sensor)?

[sjperkins]

However, if the actual data payload is msgpacked/encoded anyway, how about the following:

Here's a more concrete example:

Isn't that going to have to read over the entire file in order to apply the filter, even if you're only interested in a subset of sensors (or a single sensor)?

According to my understanding, it will only read the indexing columns of the primary key (i.e. timestamp and sensor_id) to identify the relevant rows. Note that sensor_id would probably be constructed as a categorical array.

Then, it will read those rows of the payload column out of the relevant row goups (chunks) within the parquet file.