brew install python   # installs into /opt/homebrew/bin/python3
export PATH="$(brew --prefix)/opt/python@3/libexec/bin:$PATH"   # add this to your ~/.bashrc
                                                                # so that it finds pip, python
brew install hdf5         # C-flavoured HDF5 library
brew install virtualenv
pip install --upgrade pip
pip install pip_search
pip_search tables   # first entry, version 3.8.0, "Hierarchical datasets for Python"

virtualenv ~/pytables-env
source ~/pytables-env/bin/activate
env HDF5_DIR=/opt/homebrew/opt/hdf5 pip install tables   # install PyTables
pip install netCDF4       # if you want to work with NetCDF data
pip install Pillow        # if you want to work with images
pip install wonderwords   # if you want to name stars
...
deactivate

module load netcdf    # only if you want to play with the NetCDF example below
module load hdf5/1.12.1 python/3.10.2 scipy-stack
python
>>> import tables

cd ~/tmp/pytables
/bin/rm -f *.h5 compressed.xmf
source ~/pytables-env/bin/activate
python

import tables as tt
import numpy as np

help(tt.IsDescription)          # class IsDescription inside PyTables

class Star(tt.IsDescription):   # create a child class "Star" (class inheritance)
    name     = tt.StringCol(16) # each name contains 16 8-bit numbers, could be used to store ASCII characters
    ra       = tt.Float64Col()  # Right Ascension as double-precision float
    dec      = tt.Float64Col()  # Declination as double-precision float
    white    = tt.Float32Col()  # white magnitude as single-precision float
    red      = tt.Float32Col()  # red magnitude as single-precision float
    blue     = tt.Float32Col()  # blue magnitude as single-precision float
    parallax = tt.Float32Col()  # parallax as single-precision float
    proper   = tt.Float32Col()  # proper motion as single-precision float

h5file = tt.open_file("catalogue.h5", mode="w", title="Stellar catalogue")
g1 = h5file.create_group("/", 'data', 'Observational data')    # new group under root node; title and description
table = h5file.create_table(g1, 'stars', Star, "All stars") # table stored as a node named `stars`
                                                               # while the `Star` class defines its columns
print(h5file)   # show the object tree
h5file          # show more info, including the table columns

import random
from math import asin, degrees
from wonderwords import RandomWord
r = RandomWord()
star = table.row   # get a pointer to this table's row instance
for i in range(100):
    # star['name'] = 'star' + str(random.randint(0,999))
    star['name'] = r.word(include_parts_of_speech=["nouns"], word_min_length=4, word_max_length=8)
    star['ra'] = 360 * random.random()
    star['dec'] = degrees(asin(2*random.random()-1))
    star['white'] = random.gauss(mu=0, sigma=1)
    star['red'] = random.gauss(mu=0, sigma=1)
    star['blue'] = random.gauss(mu=0, sigma=1)
    star['parallax'] = 1.e-4 * random.random()
    star['proper'] = random.random()
    star.append()   # insert this new star into the catalogue

table.shape     # nothing's been written yet
table.flush()   # flush the table's I/O buffer: write to disk, free memory resources
table.shape     # now we have data in the table
h5file.close()

import tables as tt
import numpy as np
h5file = tt.open_file("catalogue.h5", mode="a")   # append mode

table = h5file.root.data.stars
table.cols               # 8 columns from Star class
table.shape              # 10 rows
table.col('name')        # data in the "name" column
table.col('parallax')    # data in the "parallax" column
table.cols.parallax[:]   # same
for row in table.iterrows():
    print(row[:])   # print each row

min(table.col('parallax')), max(table.col('parallax'))
names = [row['name'] for row in table.iterrows() if 3e-5 <= row['parallax'] < 5e-5]
names    # star names stored as 16-element binary objects (not strings)
b'A' == b'\x41'   # true; each 8-bit number (0-255) can store a character

condition = "dec > 45"
for row in table.where(condition):
	print(row['name'], row['dec'])

condition = "(dec > 45) & (blue > 0.5)"
for row in table.where(condition):
	print(row['name'], row['dec'], row['blue'])

names = [row['name'] for row in table.where("(dec > 45) & (blue > 0.5)")]
names

sel = h5file.create_group(h5file.root, "selection", "Northern blue stars")         # add a group: name and description
h5file.create_array(sel, 'name', np.array(names), "Northern blue stars selection") # add an array of names

h5file.close()

h5dump catalogue.h5     # full content of the file as text
h5ls -rd catalogue.h5   # shorter format

import tables as tt
import numpy as np
h5file = tt.open_file("catalogue.h5", mode="a")   # append mode

for node in h5file: # root, 2 groups /data and /selection, 1 table /data/stars, 1 array /selection/name
    print(node)

for group in h5file.walk_groups():   # root and 2 groups
    print(group)

for group in h5file.walk_groups("/data"):   # 1 group under /data
    print(group)

for node in h5file.list_nodes("/data", classname='Table'):   # 1 table under /data
    print(node)
	
h5file.list_nodes("/data", classname='Table')   # more info on this table

table = h5file.root.data.stars
table.shape

table.attrs.date = "Thu, 2023-May-27 14:33"
table.attrs.temperature = 18.4
table.attrs.temp_unit = "Celsius"
table.attrs                   # list all attributes
del table.attrs.temperature   # delete an attribute
del table.attrs.temp_unit

import random
from math import asin, degrees
from wonderwords import RandomWord
r = RandomWord()
table = h5file.root.data.stars
star = table.row
for i in range(50):   # add 50 more stars to the existing table
    star['name'] = r.word(include_parts_of_speech=["nouns"], word_min_length=4, word_max_length=8)
    star['ra'] = 360 * random.random()
    star['dec'] = degrees(asin(2*random.random()-1))
    star['white'] = random.gauss(mu=0, sigma=1)
    star['red'] = random.gauss(mu=0, sigma=1)
    star['blue'] = random.gauss(mu=0, sigma=1)
    star['parallax'] = 1.e-4 * random.random()
    star['proper'] = random.random()
    star.append()   # insert this new star into the catalogue

table.shape   # 100 original stars
table.flush()
table.shape   # now 150 stars

bodies = table.cols.name
bodies[:10]   # the first 10 star names (out of 150)

bodies[:2] = [b'one', b'two']   # overwrite the first two values
bodies[:]

table.flush()   # flush the table's I/O buffer: write to disk, free memory resources
h5file.close()

wget https://bit.ly/paraviewzipp -O paraviewWorkshop.zip
unzip paraviewWorkshop.zip data/mandelbulb300.nc && /bin/rm paraviewWorkshop.zip
ls -l data/mandelbulb300.nc   # 1MB file

cd ~/tmp/pytables
source ~/pytables-env/bin/activate
python

import tables as tt
import numpy as np

import netCDF4 as nc
ds = nc.Dataset("data/mandelbulb300.nc")
print(ds)
print(ds["stability"].shape)

h5file = tt.open_file("uncompressed.h5", mode="w", title="Mandelbulb at 300^3 resolution")
g1 = h5file.create_group("/", 'g1', 'Fractal data')   # new group under root node; title and description
h5file.create_array(g1, 'stability', ds["stability"][:,:,:], "stability variable")
h5file.close()

compress = tt.Filters(complib='zlib', complevel=5)   # create the compression filter
h5file = tt.open_file("compressed.h5", mode="w", title="Mandelbulb at 300^3 resolution", filters=compress)
g1 = h5file.create_group("/", 'g1', 'Fractal data')   # new group under root node; title and description
h5file.create_carray('/g1', 'stability', obj=ds["stability"][:,:,:])   # use chunked arrays, write into g1
h5file.close()

import tables as tt
import numpy as np
compress = tt.Filters(complib='zlib', complevel=5)   # create the compression filter
h5file = tt.open_file("compressed.h5", mode="a", filters=compress)   # append mode
a = np.random.randn(100,100)   # 100^2 array from a normal (Gaussian with x0=0, sigma=1) distribution
h5file.create_carray('/g1', 'a', obj=a)   # use chunked arrays, write into g1
h5file.close()

import tables as tt
f1 = tt.open_file("compressed.h5", mode="r")   # read mode
for node in f1:   # root, 1 group g1, 2 arrays inside g1, both compressed
    print(node)

for array in f1.walk_nodes("/", "Array"):
    print('---')
    print(array)
    print(array.shape)
    # print(array[:])

f1.root.g1.a[:]                  # entire array

f1.root.g1.stability[:]          # entire array
f1.root.g1.stability[:3,:3,:3]   # 3x3x3 starting corner block
f1.root.g1.stability.shape
f1.root.g1.stability.filters
f1.root.g1.stability.title
f1.root.g1.stability.dtype
f1.root.g1.stability.size_in_memory   # uncompressed in memory
f1.close()

import tables as tt
f1 = tt.open_file("compressed.h5", mode="a")   # append mode
f1.root.g1.a.remove()
f1.close()

import tables as tt
f1 = tt.open_file("compressed.h5", mode="a")   # append mode
f1.create_soft_link("/", "stab", "/g1/stability")   # link in root called "stab" pointing to "/g1/stability"
f1.root.stab.shape
f1.root.stab[:]
f1.close()

which ptdump           # provided by PyTables
ptdump compressed.h5   # very useful command: dumps all objects from the file

<?xml version="1.0" ?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="2.0">
 <Domain>
   <Grid Name="mesh" GridType="Uniform">
     <Topology TopologyType="3DCoRectMesh" NumberOfElements="300 300 300"/>
     <Geometry GeometryType="Origin_DxDyDz">
       <DataItem Dimensions="3" NumberType="Float" Precision="4" Format="XML">
          1 2 3
       </DataItem>
       <DataItem Dimensions="3" NumberType="Float" Precision="4" Format="XML">
        0.1 0.1 0.1
       </DataItem>
     </Geometry>
     <Attribute Name="stability" AttributeType="Scalar" Center="Node">
       <DataItem Dimensions="300 300 300" NumberType="Int" Precision="4" Format="HDF">
        compressed.h5:/g1/stability
       </DataItem>
     </Attribute>
   </Grid>
 </Domain>
</Xdmf>

import tables as tt
from tables.nodes import filenode

f1 = tt.open_file('storage.h5', 'w')                         # create a new HDF5 file
fnode = filenode.new_node(f1, where='/', name='text_file')   # create a new node file inside this file
print(f1.get_node_attr('/text_file', 'NODE_TYPE'))           # it looks like a file inside our HDF5 file

fnode.write("This is a sample line.\n".encode('utf8'))
fnode.write("Write a second line.\n".encode('ascii'))
fnode.write(b"Write a third line.\n")

fnode.seek(0)   # rewind to the beginning of the file
for line in fnode:
    print(line)

fnode.close()
print(fnode.closed)   # check that it's closed
f1.close()

ls -l storage.h5    # 68K
ptdump storage.h5   # root /, file /text_file

import tables as tt
from tables.nodes import filenode
f1 = tt.open_file('storage.h5', 'r')
fnode = filenode.open_node(f1.root.text_file, 'r')
for line in fnode:
    print(line)

fnode.close()
f1.close()

import tables as tt
from tables.nodes import filenode
f1 = tt.open_file('storage.h5', 'a')                  # append mode
fnode = filenode.open_node(f1.root.text_file, 'a+')   # read and append mode
fnode.write(b"Write a fourth line.\n")
fnode.seek(0)   # rewind to the beginning of the file
for line in fnode:
    print(line)

fnode.close()
f1.close()

ls -l tuscany.avif               # 2874x2154 file
width=2874
height=2154
for num in $(seq -w 00 19); do   # crop it into twenty 300x300 random images
    echo $num
    x=$(echo "scale=8; $RANDOM / 32767 * ($width-300)" | bc)   # $RANDOM goes from 0 to 32767
    x=$(echo $x | awk '{print int($1+0.5)}')
    y=$(echo "scale=8; $RANDOM / 32767 * ($height-300)" | bc)
    y=$(echo $y | awk '{print int($1+0.5)}')
	convert tuscany.avif -crop 300x300+$x+$y small$num.png
done

from PIL import Image
image = Image.open('small00.png')
image.show()
# image.save('copy00.png')

import tables as tt
from tables.nodes import filenode
from PIL import Image
from glob import glob

f1 = tt.open_file('storage.h5', 'a')
tuscany = f1.create_group("/", 'tuscany')

for input in glob("small*.png"):
    print('copying', input)
    image = Image.open(input)
    fnode = filenode.new_node(f1, where='/tuscany', name=input.replace('.png', ''))
    image.save(fp=fnode, format='png')
    fnode.close()

f1.close()

/bin/rm -f small*png
ls -l storage.h5    # 2.2M - all our images are now stored inside this HDF5 file, along with the earlier text_file
ptdump storage.h5   # very useful command: dumps all objects from the file

import tables as tt
from tables.nodes import filenode
from PIL import Image

f1 = tt.open_file('storage.h5', 'r')

image = Image.open(fp=f1.root.tuscany.small08)
image.show()

for i in f1.root.tuscany:   # cycle through all its children
    print(i)
 
f1.root.tuscany.small00      # one possible syntax
f1.root.tuscany['small00']   # another syntax

image = Image.open(fp=f1.root.tuscany['small08'])
image.show()

import tables as tt
from tables.nodes import filenode
from PIL import Image

f1 = tt.open_file('storage.h5', 'a')
f1.root.tuscany.small19.remove()               # remove the node
f1.root.tuscany.small18.rename("last_image")   # rename the node

f1.root.tuscany             # all in one array

for i in f1.root.tuscany:   # cycle through all its children
    print(i)