Add function similar to gridboxmean including sellevel, seltime functionality to Diagnostics

[se04ber]

Proposal:
Add a gridboxmean(press,x) similar function that:

1. Averages horizontal grid cells of block size x together, similar to the functionality of CDO gridboxmean, but generalized for a data array as input (compatible with the output of IO_read_pfb() used in ParFlowDiagnostics for further functionality (for 3D and 2D).
2. Has incorporated sellevel and seltime included. (Can overgive also 3D and 4D arrays with choosen Zlayerindex, Tstepindex to make 3D or 2D).
3. Gives back heatarray (generally works with heat), but also possibility to give back/directly convert to netcdf file.

This function could be used at the end of a functionality chain in Diagnostics to average down output for different reasons f.e. achieve a smoother picture, when high res not needed and/or faster look at or visualization if desired, or to generally get a netcdf file with locally manageable size (MBit instead of GBit).

Prototype output file for netcdf=true tested seems in first cases to have similar results to cdo gridboxmean,10,10 for converted netcdf files in the order of accuracy up to e-10:

#####Code start
data = IO.read_pfb(path2file)
var= "et"
block_size=10
Zlayer = -1
Tstep = 0 #f.e. for forcing.nc files helpful
make2D = True
toNetcdf=True

#Here if part of Diagnostics class (self, ) added and nz,ny,nx can also be received through nz = self.Nz (..)

def gridboxmean_pfb(data, var, block_size=10, Zlayer=-1,Tstep=None,make2D=False,toNetcdf=True):
nz, ny, nx = data.shape
#New (x,y) dimensions
new_ny = ny // block_size + (1 if ny % block_size else 0)
new_nx = nx // block_size + (1 if nx % block_size else 0)

#Throw error if data array is 0- or 1-dimensionl
if len(data) < 2:
    raise ValueError("Data dimensions must larger than 1")
#Throw error if block_size does not fit to horizontal data shape
if nx % block_size != 0 or ny % block_size != 0:
    raise ValueError("Data dimensions must be divisible by block_size")

#Case 1: Averaging for all z-layes, give out 3D array of down-averaged horiozontal grid cells
if make2D == False:
    if len(data)==4:
        data = data[Tstep,:,:,:] # !!!!For now consider only 3D /convert to 3D if 4D

    newShape = ["nz","new_ny","new_nx"]
    reduced_data = ht.zeros((nz,new_ny,new_nx),split=-1)

    #Iterate and average over (block_size)*(block_size)d horizontal grid cells blocks 3D/for all layers
    for i in range(nz):
        for j in range(new_ny):
            y_start = j * block_size
            y_end = min((j + 1) * block_size, ny)     
            for k in range(new_nx):
                x_start = k * block_size
                x_end = min((k + 1) * block_size, nx)
            
                # Extract block and compute average
                block = data[y_start:y_end, x_start:x_end]
                block_mean = ht.mean(block)
                #write to reduced grid
                reduced_data[i,j, k] = block_mean

#Case 2: Select specific Z-layer, give out down-averaged 2D Array               
else:
    if len(data)==2:
        raise ValueError("Data dimensions already 2D")
     
    #Get specified time/iteration step
    if len(data)==4:
        data = data[Tstep,:,:,:] 
    #Get specified Z-layer
    data = data[Zlayer,:,:]
    
    newShape = ["new_ny","new_nx"]
    reduced_data = ht.zeros((new_ny,new_nx),split=-1)

    #Iterate and average over (block_size)*(block_size)d grid cell blocks 2D
    for j in range(new_ny):
            y_start = j * block_size
            y_end = min((j + 1) * block_size, ny)
            for k in range(new_nx):
                x_start = k * block_size
                x_end = min((k + 1) * block_size, nx)

                # Extract block and compute average
                block = data[y_start:y_end, x_start:x_end]
                block_mean = ht.mean(block)
                #write to reduced grid
                reduced_data[j, k] = block_mean

# Save to netcdf
if toNetcdf==True:
    ht.save_netcdf(
            data = reduced_data,
            path = str(output_file) + ".nc",
            variable= var,
            mode="w",
            dimension_names=newShape
            )
return reduced_data

####Code end