# Performance optimization processing 95,000+ NetCDF files (using xarray + dask)

I was told over at gis.stackexchange.com that this question would be a better fit here.

I already have a working code but it is very slow and my hunch is that it could be optimized. The goal here is to process 11 years of data = 11 x 12 x 30 x 24 = roughly 95,040 hourly files

On a RAID/NAS external storage I have NetCDF files from NLDAS-2 products that look like this:

<xarray.Dataset>
Dimensions:      (lon: 464, lat: 224, time: 1, bnds: 2)
Coordinates:
* lon          (lon) float32 -124.9 -124.8 -124.7 ... -67.31 -67.19 -67.06
* lat          (lat) float32 25.06 25.19 25.31 25.44 ... 52.69 52.81 52.94
* time         (time) datetime64[ns] 1992-11-09T15:00:00
spatial_ref  int32 0
Dimensions without coordinates: bnds
Data variables:
time_bnds    (time, bnds) datetime64[ns] 1992-11-09T14:00:00 1992-11-09T1...
Tair         (time, lat, lon) float32 ...
Qair         (time, lat, lon) float32 ...
PSurf        (time, lat, lon) float32 ...
Wind_E       (time, lat, lon) float32 ...
Wind_N       (time, lat, lon) float32 ...
LWdown       (time, lat, lon) float32 ...
CRainf_frac  (time, lat, lon) float32 ...
CAPE         (time, lat, lon) float32 ...
PotEvap      (time, lat, lon) float32 ...
Rainf        (time, lat, lon) float32 ...
SWdown       (time, lat, lon) float32 ...


The files are hourly and I'm aggregating them on a daily basis. Here are the important steps of my code:

I'm reading the data in chunks of one month at a time.

forcing_data = xr.open_mfdataset(file_list, chunks={'time': 24},
parallel=False, combine='nested', concat_dim='time',
autoclose=True)


I tried using parallel=True but with the files stored on a NAS I was getting reading errors different every time I re-ran the code.

I have a little validation step where I use the time_bnds variable to make sure I have the right time boundaries. I think this could be optimized when generating the files to read but I don't think it's the main problem here

idx_from, idx_to = get_time_indices_from_time_bnds(forcing_data, start_date, end_date)
forcing_data = forcing_data.isel(time=slice(idx_from, idx_to+1))


Augmenting

I need to generate new variables by performing arithmetic on the existing variables. This is pretty fast typically.

forcing_data['RH'] = compute_RH(forcing_data['PSurf'], forcing_data['Tair'], forcing_data['Qair'])
forcing_data['Tdew'] = compute_dew_point(forcing_data['RH'], forcing_data['Tair'])
forcing_data['wind_speed'] = compute_wind_speed(forcing_data['Wind_E'], forcing_data['Wind_N'])
forcing_data['wind_direction'] = compute_wind_direction(forcing_data['Wind_E'], forcing_data['Wind_N'])


Aggregating (Main Loop)

In a loop I aggregate the data month by month (I tried different time steps, it actually seemed like weekly/daily would work better in some cases) by reading the files in the way mentioned above

# get the files for each dataset type requested
agg_list = []

curr_date = start_date

while curr_date < end_date:

ds = read_forcing_dataset(files_dir, curr_date, one_month, product=product, augment=True)

# Aggregate the forcing data
agg = aggregate_forcing_data_daily(ds)

agg = agg.compute()

agg_list.append(agg)

curr_date += one_month

# Concatenate the aggregated data
agg_ds = xr.concat(agg_list, dim='time')


The aggregate_forcing_data_daily function has a similar while loop with a day increment, it basically uses agg_data = ds.isel(time=slice(idx_from, idx_to + 1)).map(agg_func, keep_attrs=True) with agg_func returning x.mean(dim='time') or x.sum(dim='time') based on the x.name. It also adds the minimum and maximum temperature for that day in the following fashion:

agg_data['Tmin'] = ds['Tair'].sel(time=slice(start_date, end_date)).min(dim='time', skipna=True, keep_attrs=True)
agg_data['Tmax'] = ds['Tair'].sel(time=slice(start_date, end_date)).max(dim='time', skipna=True, keep_attrs=True)


Performance so far:

This part of the code takes about 5 to 6 minutes per month of data on a laptop with Intel i9 @ 2.5 GHz and 32 GB of RAM, 4 processes (16 threads)

Extracting time-series for specific coordinates

for lat, lon in zip(lats,lons):
df = agg_ds.sel(lat=lat, lon=lon, method='nearest').to_dataframe()
df = function_to_rename_and_reorder_columns(df)

df.to_parquet(custom/path/to/savefile/latlon.parquet)


Performance for this step

Before I used agg = agg.compute() in the aggregation loop, this step took a lot of time but now it seems to be much faster. Prior to that the DataFrame generation step took looooot of time, many minutes if not hours. As it is it would take over a day of non-stop computation to process the 11 years dataset. Is it just the way it has to be?

I tried to play with the chunk size and the timestep for reading files and monitoring with the dask dashboard but I haven't really been able to improve further.

Like I said at the beginning this works but it feels like it could be improved. I think that the fact that the files are stored on a NAS is not the best but I cannot get around that for now. This is my first big code involving xarray and NetCDF files so I might be missing obvious things.