Zarr output too many files

[SeanBensadon]

Question

Question

Hi,

I am running an Ocean Parcels simulation with ~190,000 particles per release and tracking them for 30 days. However, I am encountering a disk quota limit on the number of files allowed when running the simulation on a supercomputer due to the zarr output format. Is there a way to work around this issue? For instance, could I export the output directly as a single NetCDF file instead of generating heaps of files in the zarr output?
I've included my code below for reference. Thank you very much for your assistance!

Supporting code/error messages

# Paste your code within this block

import pandas as pd
import pandas
import numpy as np
from datetime import timedelta
import matplotlib.pyplot as plt
import xarray as xr
import parcels
from parcels import FieldOutOfBoundError, ParticleSet, FieldSet, Field
from parcels import StatusCode
from parcels import DiffusionUniformKh

from shapely.geometry import Polygon, Point
import os
import shutil # Import shutil for removing directories

pathout = '/scratch/pawsey0106/sbensadon/OceanParcels/CoralBay/SensitivityAnalysis/'

Define release schedule with 24-hour run-up time

def create_release_times(start_time, durations, intervals):
release_times = []
run_up_time = 43.5 * 3600 # 24 hours + 19.30 in seconds
start_time += run_up_time # Add the run-up time to the initial start_time

for duration, interval in zip(durations, intervals):  # for each pair of duration and interval from the durations and intervals lists
    release_times.append(np.arange(start_time, start_time + duration * 3600, 3600))  # generate an array of times starting at start_time, ending at start_time + duration (converted to seconds), with 1-hour intervals (3600 seconds)
    start_time += duration * 3600 + interval * 3600  # update start_time to the end of the current release period plus the interval period (both in seconds)

return np.concatenate(release_times)  # combine all release time arrays in the list into a single array and return it

Loop through the years 1993-2023

for year in range(1993, 2008):
# Load dataset
ds = xr.open_dataset(f'/scratch/pawsey0106/sbensadon/OceanParcels/CoralBay/SensitivityAnalysis/FlowFM_sedtrails_gridded_{year}.nc')

# Define FieldSet
variables = {'U': 'U', 'V': 'V'}
dimensions = {'U': {'time': 'time', 'lat': 'y', 'lon': 'x'},
              'V': {'time': 'time', 'lat': 'y', 'lon': 'x'}}

fieldset = FieldSet.from_xarray_dataset(ds, variables, dimensions, mesh ='flat', interp_method={
    "U": "partialslip",
    "V": "partialslip",
})

fieldset.add_constant_field("Kh_zonal", 0.2)
fieldset.add_constant_field("Kh_meridional", 1)

# Extract longitude and latitude coordinates
lon_values = ds[dimensions['U']['lon']].values
lat_values = ds[dimensions['U']['lat']].values

# Define particle release locations
csvFile = pandas.read_csv('/scratch/pawsey0106/sbensadon/OceanParcels/CoralBay/release_locations_equal.csv')
lons = csvFile['xcoord'].values
lats = csvFile['ycoord'].values

# Land mask
landmask = fieldset.U.data[1,:,:]

landmask = (landmask == 0).astype('int') # if zeros
[xx,yy]=np.meshgrid(fieldset.U.grid.lon,fieldset.U.grid.lat)
fieldset.add_field(Field(name="landmask", data=landmask, lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat, mesh="flat", interp_method='nearest'))


# Define durations and intervals (in hours)
durations = [4, 4, 4]  # durations for release periods
intervals = [20, 20, 24*30]  # intervals for non-release periods

# Create release times
release_times = create_release_times(0, durations, intervals) # call the function create_release_times with a start_time of 0 and the durations and intervals lists

total_runtime_hours = sum(durations+intervals)  # Sum of all release periods

# Create a particle set with repeated positions and times
n_locations = len(lons) # count the number of locations based on the length of the lons array
n_times = len(release_times) # count the number of release times based on the length of the release_times array

# Create release times for each particle
particle_times = np.repeat(release_times, n_locations) # create an array by repeating the release_times array for each location
particle_lons = np.tile(lons, n_times) # create an array by repeating the lons array for each release time
particle_lats = np.tile(lats, n_times) # create an array by repeating the lats array for each release time

# Initialize ParticleSet
pset = parcels.ParticleSet(
    fieldset=fieldset,
    pclass=parcels.JITParticle,
    # pclass=parcels.ScipyParticle,
    lon=particle_lons,
    lat=particle_lats,
    time=particle_times
)

# Define output file path
output_file = pset.ParticleFile(
    name=f'/scratch/pawsey0106/sbensadon/OceanParcels/CoralBay/SensitivityAnalysis/pout/ParcelsOut_{year}.zarr',
    outputdt=timedelta(hours=1)
)


# Define Function for Kernel
def DeleteOutOfDomain(particle, fieldset, time):  
    if (particle.lon < fieldset.lon_min or
        particle.lon > fieldset.lon_max or
        particle.lat < fieldset.lat_min or
        particle.lat > fieldset.lat_max):
        particle.delete()


def DeleteInLand(particle, fieldset, time):  
    """Check whether particle is on land"""
    if (particle.lon < fieldset.lon_max and
        particle.lon > fieldset.lon_min and
        particle.lat < fieldset.lat_max and
        particle.lat > fieldset.lat_min):

        landcheck = fieldset.landmask[time, particle.depth, particle.lat, particle.lon]
        if landcheck == 1: 
            particle.delete()



# Add these as constants to the fieldset
fieldset.add_constant('lon_max', 790275)
fieldset.add_constant('lon_min', 772875)
fieldset.add_constant('lat_max', 7456000)
fieldset.add_constant('lat_min', 7429600)


kernel = pset.Kernel(parcels.AdvectionRK4) + pset.Kernel(DiffusionUniformKh) + pset.Kernel(DeleteInLand) + pset.Kernel(DeleteOutOfDomain)
# kernel = pset.Kernel(parcels.AdvectionRK4) + pset.Kernel(DiffusionUniformKh) + pset.Kernel(DeleteOutOfDomain)

pset.execute(
    kernel,
    runtime=timedelta(hours=total_runtime_hours),
    dt=timedelta(minutes=5),
    output_file=output_file
    )

ds_zarr = xr.open_dataset(f'{pathout}/pout/ParcelsOut_{year}.zarr', engine='zarr')
ds_zarr.to_netcdf(f'{pathout}/ParcelsOut_{year}.nc')

# Delete the Zarr file to free up space

shutil.rmtree(f'{pathout}/pout/ParcelsOut_{year}.zarr')

# Paste your error message within this block

[erikvansebille]

Hi @SeanBensadon, thanks for your question. You can change the size of each parr-chunk using the chunks keyword when you create your ParticleFile object, see also the "note on output chunking" at https://docs.oceanparcels.org/en/latest/examples/tutorial_parcels_structure.html#4.-Execution-and-output

Note that there was also a bug in Parcels that has impacted your workflow, where you release particles at different times. It has been patched (#1788) and will become part of the next version 3.1.1. But in any case, setting your own chunks would help

A minimum for you could be

# Define output file path
output_file = pset.ParticleFile(
    name=f'/scratch/pawsey0106/sbensadon/OceanParcels/CoralBay/SensitivityAnalysis/pout/ParcelsOut_{year}.zarr',
    outputdt=timedelta(hours=1),
    chunks=(len(pset), 1),
)

If you want, you can also set the 1 (for number of observations per chunk) in the last line to a larger number. That will speed up the output writing and result in fewer files, but may make the analysis a bit more difficult (because it could be that the last few observations are NaN)

Hope this helps!

[SeanBensadon]

Hi Erik, thanks so much for the quick response your solution has worked very well!
All the best
Sean