SWIFTSIM / swiftsimio

Python library for reading SWIFT data. Uses unyt and h5py.
GNU Lesser General Public License v3.0
15 stars 12 forks source link

Projection returns wrong shape #178

Closed Joeybraspenning closed 7 months ago

Joeybraspenning commented 8 months ago

When using the project_gas() method, sometimes the returned shape is not square. For FLAMINGO this seems to be the case for a fraction of the haloes. These do not cross the edge of the box, and the list provided to the 'region' keyword gives a square area.

Here's a working example for the 190th halo in the L1000N1800, we confirmed that for example halo 180 does not have this problem (and out of the few hundred we tried, only a handful fail). The 'flux' object returned by project_gas() has shape (64,63) in this case.

import h5py
import swiftsimio as sw
import unyt

soapfile = h5py.File('/cosma8/data/dp004/flamingo/Runs/L1000N1800/HYDRO_FIDUCIAL/SOAP/halo_properties_0077.hdf5', 'r')
mask = sw.mask("/cosma8/data/dp004/flamingo/Runs/L1000N1800/HYDRO_FIDUCIAL/snapshots/flamingo_0077/flamingo_0077.hdf5")
position = soapfile["SO/500_crit/CentreOfMass"][190] * unyt.Mpc
radius = soapfile["SO/500_crit/SORadius"][190] * unyt.Mpc
load_box = [[position[0] - radius, position[0] + radius], 
            [position[1] - radius, position[1] + radius],
            [position[2] - radius, position[2] + radius]]
mask.constrain_spatial(load_box)
halo_data = sw.load("/cosma8/data/dp004/flamingo/Runs/L1000N1800/HYDRO_FIDUCIAL/snapshots/flamingo_0077/flamingo_0077.hdf5", mask=mask)
halo_data.gas.flux = halo_data.gas.xray_luminosities.erosita_low
flux = sw.visualisation.projection.project_gas(
    halo_data, 
    resolution=64, 
    project="flux", 
    region=[position[0] - radius, position[0] + radius, position[1] - radius, position[1] + radius], 
    parallel = True
)
JBorrow commented 7 months ago

I honestly don't even know how this would happen

JBorrow commented 7 months ago

Ah - I think this must be related to the changes here where people wanted the ability to create non-square images: https://github.com/SWIFTSIM/swiftsimio/blob/16cbe91b0de4d081b71bd262baf1882c6ca48c8b/swiftsimio/visualisation/projection.py#L260.

I am guessing there is some round-off error here and we probably want to be using math.ceil instead of int which takes the floor?