TCDSolar / stixpy

STIX data analysis in python
https://stixpy.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
17 stars 20 forks source link

Incorrect orientation of maps #136

Closed hannahc243 closed 2 weeks ago

hannahc243 commented 1 month ago

I'm following the imaging demo code and I noticed that imaging maps seem to be oriented incorrectly in Helioprojective Cartesian coordinates for some test flares. I think the X,Y coordinates should be swapped and their signs flipped.

In the following example the flare should be at ~(1200, 200) in HPC from Solar Orbiter (see images from data center pipelines). The python imaging pipeline gives a location in HPC_solo = ~(-200,-1200).

Could the bug originate from the coordinates in the STIX frame? Let me know if I'm missing something.

Screenshot 2024-07-09 at 12 06 19 Screenshot 2024-07-09 at 11 54 31
samaloney commented 1 month ago

Yea could well be a mistake/bug. If your following the imaging demo script could you try changing the line below 90 -> -90 I was never sure of the sense of the rotation clockwise or counterclockwise.

header_hp = make_fitswcs_header(bp_image, hpc_ref, scale=[10, 10] * u.arcsec / u.pix,
                                rotation_angle=90 * u.deg + roll)

https://github.com/TCDSolar/stixpy/blob/ba0ef1c9d3e53946128afd2c3c117f5cbfd6fe82/examples/imaging_demo.py#L132C13-L132C32

@paolomassa I thought we verified this with a few examples?

hannahc243 commented 1 month ago

rotation by -90 deg puts it in the right quadrant, but the (X,Y) coordinates still need to be swapped.

samaloney commented 1 month ago

Ok thanks could you post the filename, time range and energy range too and I'll have a look? (maybe the updated image too?) Thanks in advance

hannahc243 commented 1 month ago

Sure!

File UID: 2211132304 time range = [2022-11-13T06:18:00, 2022-11-13T06:19:00] energy range = [6,10] keV

samaloney commented 1 month ago

Ok I've found a work use the transpose of the images e.g. Map((bp_image.T, header_hp)), the same can be achieved by swapping the u, v coords I'm not sure what happened or when - I'm wondering if it is related to IDL a[2,1] vs Python a[1,2] indexing.

max_stix.transform_to(hp_map.coordinate_frame)
<SkyCoord (Helioprojective: obstime=2022-11-13T06:18:29.252, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2022-11-13T06:18:29.252, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, m)
    (-20.27442933, 7.78303157, 9.61469376e+10)>): (Tx, Ty, distance) in (arcsec, arcsec, )
    (1194.06558915, 204.51417998, 1.)>

image

image

samaloney commented 1 month ago

So in the IDL code the important line is

visin.u # reform( xypi[0,*,*], npx2) + visin.v # reform( xypi[1,*,*], npx2)

and in python

 x[..., np.newaxis] * uv[np.newaxis, 0, :] + y[..., np.newaxis] * uv[np.newaxis, 1, :])

but I'm wondering if it should be

x[..., np.newaxis] * uv[np.newaxis, 1, :] + y[..., np.newaxis] * uv[np.newaxis, 0, :])
samaloney commented 1 month ago

Ok so I found the root cause and more importantly understand it now. Have a fix in TCDSolar/xrayvision/issues/73 reviews welcome.

hayesla commented 2 weeks ago

hey @samaloney - i've tried on the main branch and still have this issue. was the fix in #73?

hayesla commented 2 weeks ago

ah I realised fix is in xrayvision.

hayesla commented 2 weeks ago

this is now fixed and working (on main dev branch of xrayvision)