I’ve done something similar way too many times. In my case all the “bands” were usually the same geophysical parameter and I wanted all the images to share one colorbar.
You can use snappy to reproduce the steps you use to generate an image in the GUI (see attached snappy_write_image.py). You might need to process the file to extract the range of values for each band an uses those to compute a color table. This file is quite old, so will likely need some updating:
snappy_write_image.py (3.9 KB).
Recently I have started using the R terra library with NetCDF4-CF format data for similar tasks, but terra has been under development so not 100% reliable.