Exitation temperatures in case of octree clouds

[noltehenrik]

Hi, i'm using LOC for my bachelor thesis and i successfully run LOC with octree data. Right now i'm facing the issue to plot the exitation temperatures. In the documentation of the method LOC_read_Tex_3D there is a reference to the routine OT_GetCoordinatesAllV() which I cannot find anywhere in the source code. Can you provide me with this method or is there another routine which helps to plot exitation temperatures in case of octree clouds?

Thank you and best regards! Henrik

[mjuvela]

Hi,

The routines were indeed missing. I added some in the ISM repository, under the directory OT. For me this snippet works,

import sys
sys.path.append('/home/mika/GITHUB/ISM/OT/')
from OT_library import *

NX, NY, NZ, LCELLS, OFF, H, T, S, VX, VY, VZ, CHI = OT_ReadHierarchyV_LOC('OT.cloud')
x, y, z  =  OT_GetCoordinatesAllV(NX, NY, NZ, LCELLS, OFF, H, GPU=0)

and you of course need to update the path.

[noltehenrik]

Thank you, that already helped a lot! Right now i am facing the issue of actually plotting the data. My first approach was to map the coordinates from OT_GetCoordinatesAllV() onto a e.g. 512x512x512 cube and fill it with the data from the .tex file. Unfortunately the output looks like that: I also tried the routine OT_GetValueV() to get values cell by cell from my Tex vector, but the routine doesnt work for me and always throws: Here TEX_OT is the vector i get if i read my .tex-file with fromfile(). I would be very thankful if you can help me with my problem.

[noltehenrik]

P.S.: The code for my first approach:

mport sys
sys.path.append('/home/nolte/Testing/ISM-Code/ISM/OT/')
from OT_library import *
import numpy as np
from numpy import *
from matplotlib.pyplot import *

filename = 'test_CO_01-00.tex'
mp_size = 256

def create_empty_array(N):
    array = np.zeros((N, N, N), dtype=float32)
    return array

def convert_index(i):
    # map index from ot on cube with length mp_size
    return floor(i*mp_size)

NX, NY, NZ, LCELLS, OFF, H, T, S, VX, VY, VZ, CHI = OT_ReadHierarchyV_LOC('oct_test_0060.cloud')
x, y, z = OT_GetCoordinatesAllV(NX, NY, NZ, LCELLS, OFF, H, GPU=0)

x = np.apply_along_axis(convert_index, 0, x)
y = np.apply_along_axis(convert_index, 0, y)
z = np.apply_along_axis(convert_index, 0, z)

x = x.astype(int)
y = y.astype(int)
z = z.astype(int)

# sys.exit()

# read tex from octree file
fp = open(filename, 'rb')
_, _, _, _ =  np.fromfile(fp, np.int32, 4)
TEX_OT =  np.fromfile(fp, np.float32)
fp.close()

TEX_OT_sorted = create_empty_array(mp_size)
# sys.exit()

for i, value in enumerate(TEX_OT):
    # only append leaf cells
    if np.isnan(value) or value < 0.0:
        continue

    # print(value)

    x_coord = x[i]
    y_coord = y[i]
    z_coord = z[i]

    # break

    if TEX_OT_sorted[x_coord, y_coord, z_coord] == 0.0:
        TEX_OT_sorted[x_coord, y_coord, z_coord] += value
    else:
        TEX_OT_sorted[x_coord, y_coord, z_coord] = (TEX_OT_sorted[x_coord, y_coord, z_coord] + value)/2 

    # print(TEX_OT_sorted[x_coord, y_coord, z_coord])

# print(TEX_OT_sorted[256, 256, 256])

im = imshow(TEX_OT_sorted[mp_size//2, :, :])
colorbar(im, orientation = 'vertical', label = r"$T_{ex} \/ \/ [K]$")
show()

And the code which throws the error:

import sys
sys.path.append('/home/nolte/Testing/ISM-Code/ISM/OT/')
from OT_library import *
import numpy as np

filename = 'test_CO_01-00.tex'
output_map = 512

def get_normalized_cube_centers(N):
    # Generate N evenly spaced values between 0 and 1 (exclusive) for each dimension
    centers = np.linspace(1/(2*N), 1 - 1/(2*N), N)

    # Use meshgrid to create a 3D grid of center coordinates
    X, Y, Z = np.meshgrid(centers, centers, centers)

    # Stack the coordinates and reshape them into a list of centers
    normalized_centers = np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T

    return normalized_centers

NX, NY, NZ, LCELLS, OFF, H, T, S, VX, VY, VZ, CHI = OT_ReadHierarchyV_LOC('oct_test_0060.cloud')
x, y, z = OT_GetCoordinatesAllV(NX, NY, NZ, LCELLS, OFF, H, GPU=0)

TEX_OT = []

fp = open(filename, 'rb')
_, _, _, _ =  np.fromfile(fp, np.int32, 4)
TEX_OT =  np.fromfile(fp, np.float32)
fp.close()
print(TEX_OT)

coords = get_normalized_cube_centers(output_map)
TEX = np.ones((output_map, output_map, output_map), np.float32)

for x_coord, y_coord, z_coord in coords:
    val = OT_GetValueV(x_coord, y_coord, z_coord, NX, NY, NZ, LCELLS, OFF, TEX_OT)
    TEX[np.floor(x_coord*output_map)][np.floor(y_coord*output_map)][np.floor(z_coord*output_map)]

[mjuvela]

Hi,

I added one example of the plotting of octree data to my ISM repository (in the ISM/OT directory)... and made small fixes to the OT_library.py script there. The OT_GetValueV() accepts only scalar coordinate values and, being written in pure Python, would be exceedingly slow for large grids. I recommend the procedure shown in the added test_plot_tex.py script: convert coordinate vectors first to cell indices with OT_GetIndicesV() and then just index the octree arrays with the returned index vector. The routine OT_GetCoordinatesAllV() that I first suggested works the other way round, converting cell indices (in this case all of them) into coordinates. However, that is not very useful if you just want to sample the volume using a regular grid. And when you need just some cross sections for plots, you can work with coordinates in that plane, which is faster than using 3D grids or first extracting coordinates for all the cells.

[noltehenrik]

Hi,

thank you so much! Now everything works for me.