off-axis slice in spherical coordinates

[fuwentai]

I think yt is absolutely awesome and has saved me a lot of time. But I've noticed that when working with data in spherical coordinates, the existing slice functionality only offers slicing directions along 'r', 'theta', and 'phi'. However, these options do not fulfill many requirements in a spherical coordinate system. I am wondering if it's possible to implement a feature where I can define slices like off-axis in Cartesian coordinates even though the data is in spherical coordinates?

[welcome]

Hi, and welcome to yt! Thanks for opening your first issue. We have an issue template that helps us to gather relevant information to help diagnosing and fixing the issue.

[neutrinoceros]

Hi. I think there are two different questions to be answered here:

1. would such an addition be welcome in the code base ?

absolutely ! The way I see it, neither off-axis slices nor support for spherical geometries is complete while this is missing.

2. would historical contributors be willing to implement it ?

That's more difficult to say. Both features have been in the code base for years at this point and as far as I recall no one has even attempted to implement this. It'd require entirely new code so it's not something we can add incrementally from what's already there. It's also possible that it would require substantial refactors of the existing, cartesian-only, off-axis internals. It's hard to say for sure, but that's my intuition of the problem.

In short: I would absolutely love to see someone championing this feature, but that's a substantial undertaking.

[cphyc]

Coincidentally, it might be possible to achieve what you're looking for relatively easily™ thanks to #4741, but this PR would need to be fixed in the first place.

If we approximate each cell by a parallelepiped, the aforementioned PR provides machinery to compute what each one would look like from an arbitrary position. If we then iterate over all parallelepiped, we have the off-axis projection. This could work in the general case for any dataset that is tiled with parallelepiped.

[chrishavlin]

I'd also find this hugely helpful.... but don't think I have the time right now to work on it in earnest.

And for what it's worth, I've spent some time exploring ways to do this with yt as-is and depending on your needs, there are a couple approaches that can work OK:

1. embed a spherical dataset within a cartesian dataset using the Stream with callables. This essentially relies on the outer cartesian dataset to act as a pixelizer. I've been working on using this approach to handle some more general transformations with yt_xarray to avoid pre-interpolation, here's a sample notebook that defines data in geocentric coordinates and loads it within a cartesian dataset (and here's what it looks like in yt_idv).
2. use derived fields and cut regions to select data from a spherical dataset within a specified distance from an arbitrary plane (example here). This just gives you points though, not a pixelized image.

[matthewturk]

I spent some time thinking about this and I think we could do it in a reasonably straightforward way with a handful of modifications to the codebase. In essence, a cutting plane through a sphere, given a normal vector, is the same as an equatorial slice through a sphere where the sphere has been rotated such that its pole is the normal vector.

So implementing the cutting plane itself would not be terribly difficult in that sense, as we could do almost everything via simply defined offsets, which could be passed in to the routines in slice_selector.pxi. (Note that this is the slice functionality, not the cutting plane functionality.) Once this is done, we could supply these, alongwith the same offset information, to the pixelization routine.

I'm more optimistic that this is feasible than I was before reading this thread. For instance, this is what the select_cell method looks like at present:

    @cython.boundscheck(False)
    @cython.wraparound(False)
    @cython.cdivision(True)
    cdef int select_cell(self, np.float64_t pos[3], np.float64_t dds[3]) noexcept nogil:
        if self.reduced_dimensionality == 1:
            return 1
        if pos[self.axis] + 0.5*dds[self.axis] > self.coord \
           and pos[self.axis] - 0.5*dds[self.axis] - grid_eps <= self.coord:
            return 1
        return 0

Assuming we've added offset[3] info to the class itself, we could modify this to use (pos[self.axis] - self.offset[self.axis]) instead of the unmodified version, and self.coord would always be the equatorial angle, dependent on the coordinate system (i.e., 0, pi/2). Similarly, fill_mask_regular_grid could use this offset when computing the icoord value.

After writing this out, I think the best solution might be to have a rotation parameter that gets passed into slices that only means anything in the case that the geometry is spherical, and to have the slice selector query that.

[chrishavlin]

a cutting plane through a sphere, given a normal vector, is the same as an equatorial slice through a sphere where the sphere has been rotated such that its pole is the normal vector.

this seems like a promising approach! Would we be limited to the case where the sphere's origin is included in the cutting plane? I'm mostly asking out of curiosity -- I think for most applications the origin be a point in the cutting plane.

Also, would anything extra have to be done to get this to work with data selection objects? Having this modified cutting plane work with ds.region in particular would be awesome cause then you could very simply create cross-sections along great-circle minor arcs (and close #3259 ?) by using a cutting plane that includes the left and right edges of a ds.region.

[fuwentai]

thank you all guys! I am currently using an external interface provided by chrishavlin to temporarily solve this problem. First, I use ds.r[::200j,::200j,::200j] to extract a three-dimensional array in spherical coordinates, and then convert it to Cartesian coordinates for input into 'yt.' Meanwhile, I find it very interesting to implement this project in the codebase. After I obtain my Ph.D. degree, I hope I have time to make contributions in this area.