fusion-flap / flap

Fusion Library of Analysis Programs
MIT License
11 stars 5 forks source link

Slicing needs revision #73

Closed thelampire closed 3 years ago

thelampire commented 4 years ago

I am doing the following:

flap.get_data('NSTX_GPI', exp_id=139901, name='', object_name='GPI_RAW')

flap.list_data_objects('GPI_RAW') returns the following:


GPI_RAW(data_source:"NSTX_GPI" exp_id:"139901") data_title:"NSTX GPI data" shape:[59650,64,80][no error] Data name:"Signal", unit:"Digit" Coords: Time s [<R. symm.>] Start: 2.500E-01, Steps: 2.515E-06 Sample n.a. [<R. symm.>] Start: 0.000E+00, Steps: 1.000E+00 Image x Pixel [<R. symm.>] Start: 0.000E+00, Steps: 1.000E+00 Image y Pixel [<R. symm.>] Start: 0.000E+00, Steps: 1.000E+00 Device R m [<R. symm.>] Start: 1.403E+00, Steps: 3.718E-03, -7.782E-04 Device z m [<R. symm.>] Start: 7.054E-02, Steps: 1.809E-04, 3.066E-03

I want to slice the data along Device R and Device z coordinates, cut out a ROI in spatial coordinates. According to the specification, this should be working, but it fails when I am doing the following:

flap.slice_data('GPI_RAW', slicing={'Device z':flap.Intervals(0.15,0.2)})

It returns with the following error message: raise ValueError("No data in slicing range.")

Obviously there is data in the range, the vertical range is from 0.070544312m to 0.32413751674m There is an area where there is not data in that range, but that shouldn't matter when slicing. Also, if I choose a spatially rectangular region like the following: flap.slice_data('GPI_RAW', slicing={'Device z':flap.Intervals(0.15,0.2), 'Device R':flap.Intervals(1.45,1.5)})

It fails with the same error message.

If I plot a sliced "signal" e.g. the following: flap.plot('GPI_RAW', plot_type='xy', slicing={'Device z':0.2, 'Time':flap.Intervals(0.3245,0.3255),'Image x':12}, axes=['Time']) it doesn't fail, but returns an obviously badly sliced plot.

The area of the measurement looks like this:

image

Could you tell me if I am doing something wrong? I really need working slicing for my research.

thelampire commented 4 years ago

Also, I would like to give another example. I noticed this error when i was plotting poloidally neighboring channels in pixel coordinates. Then I wanted to find out the velocity of a certain filament by plotting 5-6 spatially sliced channels. Let me illustrate:

Slicing along image x, image y coordinates: for i in range(8): flap.plot('GPI', plot_type='xy', slicing={'Image y':i*10 , 'Time':flap.Intervals(0.3245,0.3255), 'Image x':flap.Intervals(0,24)}, summing={'Image x': 'Mean'})

Figure_2

Slicing along device z and image y coordinates: for i in range(10): flap.plot('GPI', plot_type='xy', slicing={'Device z':(i*0.02+0.1), 'Time':flap.Intervals(0.3245,0.3255), 'Image x':flap.Intervals(0,24)}, summing={'Image x': 'Mean'}) Figure_4

If one takes a look at the filament at 0.3246s then there is an obvious image y direction movement, but on the spatially sliced image, it is not there at the left. Some channels seem OK, but the bottom green and blue ones are not ( there is signal there.

Comments:

thelampire commented 4 years ago

The problem is the following: Coordinates are equidistant, but they don't form a rectilinear grid. They are linearly dependent on image x and image y coordinates. The code assumes that simple slicing is possible, but it is not in the way it is implemented.

thelampire commented 4 years ago

Also, Flux r interpolation fails due to some bogus non-monotonicity error. I did calculate that coordinate from the EFIT and I am doing interpolation both in spatial and temporal coordinates successfully. I think slicing needs a major revision, there are just too many available, but untested options. I don't think I am doing anything too crazy. Anyways, here is a data object to play with: https://drive.google.com/file/d/11Uzxe6Jm4noUxBcZabKZKY0hE4YWqsdi/view?usp=sharing

thelampire commented 4 years ago

If I add the device coordinates as a matrix and not use the equidistant option with the steps, the slicing is working. However, the 1D time tract plotting still returns incorrect time series. It doesn't matter whether I am doing the slicing along the equidistant coordinates or the non-equidistant ones.

I added the new coordinates like this:

r_coordinates=np.zeros([64,80])
z_coordinates=np.zeros([64,80])
for i_x in range(64):
    for i_y in range(80):
        r_coordinates[i_x,i_y]=coeff_r[0]*i_x+coeff_r[1]*i_y+coeff_r[2]
        z_coordinates[i_x,i_y]=coeff_z[0]*i_x+coeff_z[1]*i_y+coeff_z[2]

coord.append(copy.deepcopy(flap.Coordinate(name='Device R noneq.',
                                        unit='m',
                                        mode=flap.CoordinateMode(equidistant=False),
                                        values=r_coordinates,
                                        shape=r_coordinates.shape,
                                        dimension_list=[1,2]
                                        )))

coord.append(copy.deepcopy(flap.Coordinate(name='Device z noneq.',
                                       unit='m',
                                       mode=flap.CoordinateMode(equidistant=False),
                                       values=z_coordinates,
                                       shape=z_coordinates.shape,
                                       dimension_list=[1,2]
                                       )))
thelampire commented 4 years ago

I tried the following:

The purple spiked one is the correct signal, the others are offset by at most 0.6mm, one pixel is around 3mm.

The following plot shows that the neighboring pixels have a very similar signal:

Screen Shot 2019-12-15 at 6 29 56 PM

Does anyone know what's going on?

thelampire commented 4 years ago

When I tried the interpolation with the current coordinates, the non-monotonic error came up again even though everything is linear.

thelampire commented 4 years ago

As far as I can see it, the interpolation is not even done in 2D but linearly first in the first coordinate then the second. Should I use multislice somehow? As far as I can see in Sandor's email, I shouldn't be using multislice, this should be working. How do we fix this?

thelampire commented 3 years ago

Closing the issue as it hasn't come up again and it was only a misinterpretation.