Load scientific array data (NetCDF and GRIB2)

Applies to VAST Cluster 5.5. Examples use the VAST DataBase Python SDK 2.1.0.

NetCDF and GRIB2 files hold multi-dimensional arrays, such as temperature by time, latitude and longitude. A VAST DataBase table holds rows. This guide turns one into the other with xarray: each grid point becomes a row, and each dimension becomes a column you can filter on. It loads a NetCDF file one time step at a time so memory stays flat, loads in parallel, queries a region, and then loads a real GRIB2 weather field.

Run the code blocks in order in one Python session. Below each block is the output it produced on a VAST Cluster 5.5 system.

Before you start

You need the connection settings from the quickstart in your environment, and these packages:

pip install vastdb xarray netCDF4 cfgrib eccodes

netCDF4 reads NetCDF files. cfgrib reads GRIB files, and the eccodes package includes the ecCodes library it needs.

Connect, and create the schema for this guide:

import concurrent.futures
import os
import tempfile
import threading
import urllib.request

import numpy as np
import pandas as pd
import pyarrow as pa
import pyarrow.compute as pc
import vastdb
import xarray as xr

session = vastdb.connect()
BUCKET = os.environ["VASTDB_BUCKET"]
SCHEMA = "guide_scientific"
workdir = tempfile.mkdtemp()

with session.transaction() as tx:  # start from an empty schema, so you can run the guide again
    old = tx.bucket(BUCKET).schema(SCHEMA, fail_if_missing=False)
    if old:
        for table in old.tables():
            table.drop()
        old.drop()
with session.transaction() as tx:
    tx.bucket(BUCKET).create_schema(SCHEMA)
print("Schema ready:", SCHEMA)

Schema ready: guide_scientific

NetCDF

Make a small NetCDF file to load: 2-metre temperature every hour for a day, on a 2.5-degree global grid. Your own files will be bigger, but they have the same shape: one or more variables over named dimensions.

times = pd.date_range("2026-09-16", periods=24, freq="h")
lats = np.arange(-90, 90.1, 2.5)
lons = np.arange(0, 360, 2.5)
hour = np.arange(24)[:, None, None]
by_lat_and_hour = (288 - 30 * np.abs(np.sin(np.radians(lats)))[None, :, None]
                   + 5 * np.sin(2 * np.pi * hour / 24))
t2m = np.broadcast_to(by_lat_and_hour, (len(times), len(lats), len(lons))).astype("float32")

ds = xr.Dataset({"t2m": (("time", "lat", "lon"), t2m, {"units": "K"})},
                coords={"time": times, "lat": lats, "lon": lons})
nc_path = os.path.join(workdir, "t2m.nc")
ds.to_netcdf(nc_path)
print(dict(ds.sizes), "->", ds.sizes["time"] * ds.sizes["lat"] * ds.sizes["lon"], "grid points")

{'time': 24, 'lat': 73, 'lon': 144} -> 252288 grid points

From an array to rows

Open the file. xarray reads it lazily, so nothing is loaded into memory until you ask for it:

ds = xr.open_dataset(nc_path)
print(list(ds.data_vars), {name: str(ds[name].dtype) for name in ds.coords})

['t2m'] {'time': 'datetime64[ns]', 'lat': 'float64', 'lon': 'float64'}

Each grid point becomes a row, with one column per dimension and one per variable. Convert one time step at a time, so a big file never has to fit in memory:

grid_schema = pa.schema([
    ("time", pa.timestamp("s")),
    ("lat", pa.float32()),
    ("lon", pa.float32()),
    ("t2m", pa.float32()),
])

def time_step_rows(ds, i):
    frame = ds.isel(time=i).to_dataframe().reset_index()
    return pa.Table.from_pandas(frame[grid_schema.names], preserve_index=False).cast(grid_schema)

first = time_step_rows(ds, 0)
print(first.num_rows, "rows per time step")
print(first.slice(0, 2).to_pylist())

10512 rows per time step
[{'time': datetime.datetime(2026, 9, 16, 0, 0), 'lat': -90.0, 'lon': 0.0, 't2m': 258.0}, {'time': datetime.datetime(2026, 9, 16, 0, 0), 'lat': -90.0, 'lon': 2.5, 't2m': 258.0}]

The times come out of xarray in nanoseconds, so cast() stores them in seconds: smaller, and precise enough for hourly data. Store coordinates as plain numeric columns. Then a filter on a region or a time window is pushed down to the cluster like any other predicate; see querying.

Load it

Create the table and insert each time step:

with session.transaction() as tx:
    tx.bucket(BUCKET).schema(SCHEMA).create_table("t2m", grid_schema)

loaded = 0
for i in range(ds.sizes["time"]):
    batch = time_step_rows(ds, i)
    with session.transaction() as tx:
        loaded += len(tx.bucket(BUCKET).schema(SCHEMA).table("t2m").insert(batch))
print("loaded", loaded, "rows")

loaded 252288 rows

insert() returns the row IDs it assigned, so counting them costs nothing. To check a big load afterwards, count with SQL through the VAST Query Engine, SELECT count(*) FROM "<bucket>/<schema>".t2m, which returns only the number. Reading a column to count it sends every row to your client, and get_stats() can lag behind recent writes; see write path.

For real data volumes, choose the table structure before you load. Most queries on gridded data filter on time and location, so a sorted table with time, lat and lon as sort columns, or a partitioned table by day of time, lets those queries skip data. Sorting happens in blocks of 512K rows, so a table this small isn't sorted at all; see write path.

Load in parallel

Loading many time steps, or many files, goes faster with several inserts running at once. The NetCDF and HDF5 libraries aren't safe to call from several threads at the same time, so read the file in one thread, and run only the inserts in parallel. Give each worker thread its own session:

local = threading.local()

def insert_batch(batch):
    if not hasattr(local, "session"):
        local.session = vastdb.connect()
    with local.session.transaction() as tx:
        return len(tx.bucket(BUCKET).schema(SCHEMA).table("t2m_parallel").insert(batch))

with session.transaction() as tx:
    tx.bucket(BUCKET).schema(SCHEMA).create_table("t2m_parallel", grid_schema)

loaded, pending = 0, set()
with concurrent.futures.ThreadPoolExecutor(max_workers=4) as pool:
    for i in range(ds.sizes["time"]):
        pending.add(pool.submit(insert_batch, time_step_rows(ds, i)))  # read here, insert in a worker
        if len(pending) >= 8:  # keep at most 8 batches in memory
            done, pending = concurrent.futures.wait(pending, return_when=concurrent.futures.FIRST_COMPLETED)
            loaded += sum(f.result() for f in done)
    loaded += sum(f.result() for f in pending)
print("inserted", loaded, "rows")

inserted 252288 rows

The limit on pending batches keeps memory flat however many time steps there are. When converting the arrays is the slow part, use processes instead, one file per process: put the worker function in a module and run it with concurrent.futures.ProcessPoolExecutor. Each process then opens its own file and its own session, so nothing is shared.

Query a region

Find the mean temperature over the tropics at noon. The time and latitude filter runs on the cluster, so only the matching rows come back:

with session.transaction() as tx:
    t = tx.bucket(BUCKET).schema(SCHEMA).table("t2m")
    noon = pd.Timestamp("2026-09-16 12:00")
    rows = t.select(columns=["lat", "t2m"],
                    predicate=(t["time"] == noon) & t["lat"].between(-20, 20)).read_all()
print(rows.num_rows, "rows, mean", round(pc.mean(rows["t2m"]).as_py() - 273.15, 2), "°C")

2448 rows, mean 9.37 °C

GRIB2

GRIB2 is the format most weather models publish in. A GRIB2 file is a sequence of messages, one field each, such as temperature at 2 metres for one time. Weather services usually publish an index file (.idx) next to each GRIB2 file, which gives every message's byte offset, so you can download only the fields you need.

Fetch one field

Download the 2-metre temperature from one NOAA Global Forecast System (GFS) analysis, from the NOAA open data bucket on AWS. It reads the index, finds the field's byte range, and fetches only that range:

gfs = ("https://noaa-gfs-bdp-pds.s3.amazonaws.com/"
       "gfs.20241211/00/atmos/gfs.t00z.pgrb2.1p00.f000")
index = urllib.request.urlopen(gfs + ".idx", timeout=60).read().decode().splitlines()

for n, line in enumerate(index):
    if ":TMP:2 m above ground:" in line:
        start = int(line.split(":")[1])
        end = int(index[n + 1].split(":")[1]) - 1
        print(line)
        break

request = urllib.request.Request(gfs, headers={"Range": f"bytes={start}-{end}"})
grib_path = os.path.join(workdir, "gfs_t2m.grib2")
with open(grib_path, "wb") as f:
    f.write(urllib.request.urlopen(request, timeout=120).read())
print(os.path.getsize(grib_path), "bytes downloaded")

580:34661063:d=2024121100:TMP:2 m above ground:anl:
74048 bytes downloaded

Open it and fix the types

Open the file with the cfgrib engine. indexpath="" stops cfgrib from writing an index file next to it:

grib = xr.open_dataset(grib_path, engine="cfgrib", backend_kwargs={"indexpath": ""})
frame = grib.to_dataframe().reset_index()
print(frame.dtypes.to_dict())
print(len(frame), "rows")

{'latitude': dtype('float64'), 'longitude': dtype('float64'), 'time': dtype('<M8[ns]'), 'step': dtype('<m8[ns]'), 'heightAboveGround': dtype('float64'), 'valid_time': dtype('<M8[ns]'), 't2m': dtype('float32')}
65160 rows

step (the forecast lead time) is a duration, and the cluster doesn't accept duration columns (data types). Store it as whole seconds in an integer column, so it stays a number you can filter and sort on. The coordinates that are the same for every row of a field (time, step, heightAboveGround and valid_time) are still worth keeping as columns: once you load more fields and forecast times, they tell the rows apart.

frame["step_s"] = (frame.pop("step") // pd.Timedelta(seconds=1)).astype("int64")
gfs_schema = pa.schema([
    ("time", pa.timestamp("s")),
    ("valid_time", pa.timestamp("s")),
    ("step_s", pa.int64()),
    ("heightAboveGround", pa.float32()),
    ("latitude", pa.float32()),
    ("longitude", pa.float32()),
    ("t2m", pa.float32()),
])
gfs_rows = pa.Table.from_pandas(frame[gfs_schema.names], preserve_index=False).cast(gfs_schema)

with session.transaction() as tx:
    table = tx.bucket(BUCKET).schema(SCHEMA).create_table("gfs_t2m", gfs_schema, fail_if_exists=False)
    table.insert(gfs_rows)
print("loaded", gfs_rows.num_rows, "rows")

loaded 65160 rows

Query it

Find the coldest and warmest points in the analysis, in degrees Celsius:

with session.transaction() as tx:
    t = tx.bucket(BUCKET).schema(SCHEMA).table("gfs_t2m")
    rows = t.select(columns=["latitude", "longitude", "t2m"]).read_all()

for label, i in (("coldest", pc.index(rows["t2m"], pc.min(rows["t2m"])).as_py()),
                 ("warmest", pc.index(rows["t2m"], pc.max(rows["t2m"])).as_py())):
    r = rows.slice(i, 1).to_pylist()[0]
    print(f"{label}: {r['t2m'] - 273.15:.1f} °C at lat {r['latitude']}, lon {r['longitude']}")

coldest: -46.9 °C at lat 72.0, lon 324.0
warmest: 35.4 °C at lat -25.0, lon 137.0

To load a whole GRIB2 file with many fields, open each kind of field separately: cfgrib needs backend_kwargs={"filter_by_keys": {...}} (for example {"typeOfLevel": "surface"}) when a file mixes fields on different levels.

Clean up

grib.close()
ds.close()
with session.transaction() as tx:
    schema = tx.bucket(BUCKET).schema(SCHEMA)
    for table in schema.tables():
        table.drop()
    schema.drop()
print("Dropped", SCHEMA)

Dropped guide_scientific