I was wondering whether there’s any plan to release a storage-conscious version of the catalog files — optimized to use as little disk space as possible while still providing maximum utility to the community.
Experimenting locally, I’ve been able to repack the DIA and object light curves into HDF5 at roughly 600 GB versus the current ~5.1 TB — about an 8× reduction — with no loss of photometric precision and only the removal of exactly-derivable or rarely-used columns. Breaking that down: the DIA collection goes from ~1.2 TB to ~150 GB, and the object collection from ~3.9 TB to ~480 GB. A few things that stand out:
Redundant photometry. psfMag, psfMagErr, and psfMagErr_corrected are exact functions of the corresponding fluxes — I verified psfMag = −2.5·log₁₀(psfFlux) + 31.4 reproduces the stored values to <2×10⁻⁶ mag across several million sources. Those three columns alone are ~25% of every forced-source file, about 600 GB across both collections, storing something any user can recompute in one line.
Per-epoch storage of per-visit quantities. midpointMjdTai and band are properties of the visit, not of the object-epoch — in one partition I checked, 568 unique visits were carrying 1.7 million copies of the same timestamps. Splitting into a small visit/exposure table plus a visit_idx per epoch removes that duplication entirely, and has the nice side effect of making per-object access a contiguous slice rather than a full column-chunk decode.
Compression settings. The files are SNAPPY-compressed, which achieves a ratio of about 1.01 on float32 photometry — essentially nothing. Switching to zstd with BYTE_STREAM_SPLIT encoding on the float columns gives 30–37% for free, with no schema change and slightly faster reads.
Error columns. psfFluxErr is float32, but float16 reproduces it to a median relative error of 1.7×10⁻⁴ with no overflow — nobody needs an error bar to four decimal places.
Individually these are small; together they’re most of the volume. It’s no big deal at a few TB, but I can imagine it paying big dividends as the survey grows — both in storage and in download times for users pulling light curves at scale. I’d frame it as a light-curve-optimized companion product alongside the full catalogs rather than a replacement, since the compact form necessarily drops columns some users will want.
Good question — I should have been specific. The comparison is against the HATS/LSDB parquet distribution of DP2, i.e. the bulk light-curve files, not TAP results or the Butler. Concretely:
dia_object_collection/dia_object_lc — 8,393 partitions, 232M objects, 1,197 GB on disk
object_collection/object_lc — 8,840 partitions, 786M objects, 3,911 GB on disk
Both are Parquet v2.6, SNAPPY-compressed, one row group per partition file, with the light curves in a nested struct-of-lists column (diaObjectForcedSource / objectForcedSource). My numbers are du on those trees versus du on the HDF5 I generate from them.
One important caveat so the comparison isn’t misread as apples-to-apples: my repack keeps the forced-source light curves (time, band, flux, flux error, quality flags, visit/detector provenance) and drops the rest — diaSource (96 columns of difference-image measurements), the difference-image fluxes, and the per-object scalar catalog, which is 67% of the object collection’s volume on its own. So a fair statement of the claim is:
Of the reduction, the part that is strictly lossless is the derived magnitude columns (exactly recomputable from flux) plus better compression settings — roughly 45% of the DIA tree.
The remainder comes from dropping columns my use case doesn’t need, which is a projection rather than a compression.
That’s why I’d pitch it as a light-curve-optimized companion product rather than a replacement. I retain visit_id per epoch precisely so anything dropped can be joined back against the full catalogs.
Happy to share the per-column byte accounting if useful — it’s just Parquet footer metadata, so it’s cheap to reproduce on any partition.
When it comes to “as the survey grows” you should be aware that the intent in the future is to provide per-source TAIs, taking the shutter motion into account.
It is possible that the project could look at providing them as a base time per visit and a delta per source, to make the delta more compact and compressible. That is still future work.
Good to know, thanks! A per-source delta (or even an analytic function users could evaluate themselves) would be excellent.
That’s essentially the approach I took for the barycentric correction: rather than a per-source value, I compute one BJD per (visit, detector), evaluated at the centroid of the sources landing on that detector. The residual error is then bounded by the detector’s angular size — I measure a worst case of ~1.2 s against a per-source astropy calculation.
For scale, the barycentric correction varies by ~8.7 s/degree, so across the 3.5° field of view a single field-centre value is off by up to ~15 s at the edge. The part that actually matters for time-series work is that the position-dependent term has an annual modulation, so over a ~240 d baseline it drifts by ~25 s peak-to-peak — nearly half a cycle for a 1-minute period, though negligible above an hour. Per-detector brings that to ~1 s, and per-source removes it entirely.
One practical note if the project does pursue a compact delta: the delta needs to be defined against a low-magnitude epoch rather than in the full JD frame. A float32 storing a JD-frame quantity (~2.46×10⁶) has a ULP of ~340 s, but a float32 storing just the correction (bounded by ±500 s) is good to ~30 µs — and since float64 at JD magnitude only carries ~47 µs anyway, a float32 delta plus a float64 base time is lossless in practice. That combination is very compact and very compressible.
In any case, I think a great way to boost community use of the data is getting creative about squeezing the essentials into as little disk space as possible — especially with storage costs where they are. I’d imagine it would take real pressure off the network infrastructure too.
The natural base time would of course be something like the midpoint TAI, so your suggestion was implicit there, hence the “more compact”. It would probably even be acceptable to use a short in milliseconds. I think the project would rather not use a float16 as the original authoritative source, as it’s not universally supported by standards.
Yes, good point! And agreed — I can see the project not wanting float16 as the authoritative source, though a lighter derived release could certainly use it.
On that note, another thought for saving a lot of space: rather than storing errors in flux units, store them as a fractional error that you multiply the flux by. That decouples the error’s dynamic range from the flux scale, which is what forces the wide float type in the first place.
My thoughts (from HATS/LSDB team), roughly in the order of importance:
Snappy vs ZSTD - that is an easy win, I actually thought we were using ZSTD (pull here), but it seems that it somehow was not applied in our pipeline. We are investigating and this is a valuable feedback!
Extra columns for magnitudes - we actually ADD these intentionally, exactly so that an user does not have to do a conversion by themselves. The idea is that an user, if they dont want those, they can just not load them because Parquet enables column reads so easily. So, in order to create a custom smaller catalog one could run write_catalog without unwanted columns to write down smaller catalog aimed at your purpose.
2.a) I do take the point that while Observatory does not have a way to download bulk data, there is no easy way to do write_catalog option that I described above.
Same as above for midpointMjdTai. We intentionally do not want a normal user to have to be joining tables which is somewhat complex operation, and not very performant. Again, if you want, you can remove midpointMjdTai , and use visit_table directly, which we provide in HATS directory too.
I agree about points on deltas, that is a neat idea. Just to be clear, midpointMjdTai in HATS catalog is float32, which has lower precision, it does round up to few minutes. If you need precise times, go via butler or pull original midpointMjdTai in float 64.
HDF5 has its own set of pros and cons, and we are not actively thinking about moving away from Parquet.
P.S. The original post would be easier to read for humans if the AI input was more heavily edited.
Thanks, this is helpful! And sorry about the un-edited AI output–it’s been working a lot with me on this and I’ve been going back and forth with it so much that perhaps I found it easy to read as I was already in the weeds with it
Is there a way to run write_catalog on the Rubin side of things before pulling data, or could that be set up with an accompanying tutorial so that people aren’t pulling over mostly bits that they don’t need? It’s easy for me to compress locally once pulled over, and wasn’t much of an issue this time around since it was only 5 TB of data, but I’m more concerned when it grows in scale. The stuff I work on by necessity runs over every single lightcurve, so trying to be sure that things are manageable with coming data releases.
I completely agree that what we did for DP1 and DP2 is not going to work for what we expect for DR1. It is a subject of internal discussions.
Is there a way to run write_catalog on the Rubin side of things before pulling data, or could that be set up with an accompanying tutorial so that people aren’t pulling over mostly bits that they don’t need?
No, I dont think that is possible to do it easily on RSP at the moment. I will also see if we can enable something even before DR1
Great to hear that you’re already thinking about this internally! Would be super enthusiastic to be involved in any way I can–I’m very invested in being able to pull the data over locally in a manageable way. Btw, any updated estimates on the final catalog sizes based on the depths/density of objects seen in DP2 would also be super useful I think for many of us (e.g. in knowing how more storage we need to plan for–my own naive scaling is that with the lean hdf5 storage format I’m using and extrapolating from the densest DP2 fields to the rest of the sky and 10 years worth of epochs, it’d be around a third of a petabyte to store the lightcurves after 10 years, which is actually not so scary sounding )
Btw, some really nice stuff popped out of my search of DP2–can’t wait to see how this scales to the full sky. (part of why I’m so enthusiastic about planning for DR1 and beyond)
Btw, any updated estimates on the final catalog sizes based on the depths/density of objects seen in DP2 would also be super useful I think for many of us
Object lightcurve, all together, around 1.2 PB. Under some relatively light restrictions (compressions, removing some columns; option A in the ticket) easy to drop to 550 TB. I think 300 to 400 TB at the end for just object lightcurves, with some choices being made, is very realistic.
I just wanted to add some points on top of Neven’s as an engineer on the LSDB team.
Thanks for pointing out the compression type, that was a bug that slipped through in our pipeline as we tried to get our DP2 catalogs out for the release date.
In terms of the per visit quantities, in parquet we take advantage of dictionary and run length encoding, so on disk those repeated values should have a minimal impact on the file size. Also with parquet storing data as columnar, a single object’s data is always spread across separate column chunks rather than as a contiguous byte range, so that particular per-object access optimization doesn’t really translate for us. The benefit we get from that columnar access is that for a use case like yours where you don’t need the derived columns, you can download or read only the columns you want without needing a separate optimized release. Although in this particular case where the extra columns are within our nested lightcurves, pyarrow currently doesn’t support partial nested column reading, but we are working on that!
Using byte stream split for float columns without repeated values such as fluxes is a good suggestion though, thank you! We will be looking into using that. I think between ZSTD and better optimizing dictionary vs BSS encoding per column we can get most of the benefits you’ve outlined.
