from netCDF4 import Dataset
import numpy as np

# create netCDF file
rootgrp = Dataset('geo2.nc', 'w', format="NETCDF3_CLASSIC")

dSize = 30
timeSize = 15

# 3 spatial dimensions, 4th dimension is time
xd = rootgrp.createDimension("x", dSize)
yd = rootgrp.createDimension("y", dSize)
zd = rootgrp.createDimension("z", dSize)
timed = rootgrp.createDimension("time", timeSize)

times = rootgrp.createVariable("time", np.float32, ("time",))
times.units = "seconds"		# Paraview requires units attribute to use the 4th dimension as time
times[:] = range(timeSize)

# 4D variable, initialize to all 0
xyzGrid = rootgrp.createVariable('xyzGrid', np.float32, ('time','z','y','x'), fill_value=0.)

xstart = 0
ystart = 0
zstart = 0
cubeSize = dSize/timeSize

# move a cube of data from 1 corner of the space to the diagonally opposite corner
# In Paraview, open the netCDF file and choose NetCDF files generic and CF conventions
# Volume rendering in Paraview is effective

for t in range(len(timed)):
    print(t)
    for k in range(zstart, zstart + cubeSize):
        for j in range(ystart, ystart + cubeSize):
            for i in range(xstart, xstart + cubeSize):
                xyzGrid[t, k, j, i] = float(max( (t + k + j + i) , 6))

    xstart += cubeSize
    ystart += cubeSize
    zstart += cubeSize

rootgrp.close()