--- title: "Load scientific array data (NetCDF and GRIB2)" slug: "dev-database-guides-load-scientific-data" description: "Transform NetCDF and GRIB2 files into structured data with this guide, optimizing memory and querying for weather analysis efficiently." updated: 2026-10-30T18:19:11Z published: 2026-10-30T18:19:11Z canonical: "docs.vastdata.com/dev-database-guides-load-scientific-data" --- > ## Documentation Index > Fetch the complete documentation index at: https://docs.vastdata.com/llms.txt > Use this file to discover all available pages before exploring further. # 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](https://docs.xarray.dev): 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](/develop/docs/dev-database-quickstart#2-connect) 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](/develop/docs/dev-database-guides-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](/develop/docs/dev-database-integrations-adbc), `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](/develop/docs/dev-database-concepts-architecture-write-path#table-statistics). 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](/develop/docs/dev-database-guides-data-modeling#sorted-tables) with `time`, `lat` and `lon` as sort columns, or a [partitioned table](/develop/docs/dev-database-guides-data-modeling#partitioned-tables) 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](/develop/docs/dev-database-concepts-architecture-write-path#sorted-tables). ### 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('