Xarray-Reverse-Interpolation (auf Koordinaten, nicht auf Daten)

8

Ich habe ein folgendes DataArray

arr = xr.DataArray([[0.33, 0.25],[0.55, 0.60],[0.85, 0.71],[0.92,0.85],[1.50,0.96],[2.5,1.1]],[('x',[0.25,0.5,0.75,1.0,1.25,1.5]),('y',[1,2])])

Dies ergibt die folgende Ausgabe

<xarray.DataArray (x: 6, y: 2)>
array([[0.33, 0.25],
       [0.55, 0.6 ],
       [0.85, 0.71],
       [0.92, 0.85],
       [1.5 , 0.96],
       [2.5 , 1.1 ]])
Coordinates:
  * x        (x) float64 0.25 0.5 0.75 1.0 1.25 1.5
  * y        (y) int32 1 2

oder der Einfachheit halber unten mit x und Ausgabe (z) nebeneinander sortiert.

x         z (y=1)   z(y=2)
0.25      0.33      0.25
0.50      0.55      0.60
0.75      0.85      0.71
1.00      0.92      0.85
1.25      1.50      0.96
1.50      2.50      1.10

Die Daten, die ich habe, sind das Ergebnis mehrerer Eingabewerte. Einer davon ist der x-Wert. Es gibt mehrere andere Dimensionen (z. B. y) für andere Eingabewerte. Ich möchte wissen, wann mein Ausgabewert (z) größer als 1,00 wird, wobei die anderen Dimensionen festgehalten werden und der x-Wert variiert wird. Im obigen zweidimensionalen Beispiel möchte ich die Antwort [1.03 1.32] erhalten. Weil ein Wert von 1,03 für x 1,00 für z ergibt, wenn y = 1 ist, und ein Wert von 1,32 für x 1,00 für z ergibt, wenn y = 2 ist.

Bearbeiten: Da die Ausgabe z mit zunehmendem x wächst, gibt es nur einen Punkt, an dem z 1.0 als Ausgabe hat.

Gibt es eine effiziente Möglichkeit, dies mit xarray zu erreichen? Meine eigentliche Tabelle ist viel größer und hat 4 Eingänge (Dimensionen).

Vielen Dank für jede Hilfe!

Hoogendijk
quelle

Antworten:

4

xarray hat hierfür eine sehr praktische Funktion: xr.interpSie führt eine stückweise lineare Interpolation eines xarrays durch.

In Ihrem Fall können Sie damit eine stückweise Interpolation der Punkte (x, y1) und (x, y1) erhalten. Sobald dies erledigt ist, müssen Sie nur noch den Wert Ihres interpolierten xArrays, der dem Abschlusswert Ihres interpolierten y1/y2/..Arrays zugeordnet ist, auf die Zielnummer (in Ihrem Beispiel 1,00) abrufen.

So könnte das aussehen:

y_dims = [0, 1,] 
target_value = 1.0
# create a 'high resolution` version of your data array:
arr_itp = arr.interp(x=np.linspace(arr.x.min(), arr.x.max(), 10000))
for y in y_dims:
    # get the index of closest data
    x_closest = np.abs(arr_itp.isel(y=y) - target_value).argmin()
    print(arr_itp.isel(y=y, x=x_closest))

>>> <xarray.DataArray ()>
>>> array(0.99993199)
>>> Coordinates:
>>>     y        int64 1
>>>     x        float64 1.034
>>> <xarray.DataArray ()>
>>> array(1.00003)
>>> Coordinates:
>>>     y        int64 2
>>>     x        float64 1.321


Dies funktioniert zwar, ist jedoch kein wirklich effizienter Weg, um das Problem anzugehen, und hier sind zwei Gründe, warum nicht:

  1. Mit xr.interp wird das gesamte DataArray stückweise interpoliert. Wir brauchen jedoch immer nur die Interpolation zwischen den beiden Punkten, die Ihrem Zielwert am nächsten liegen.
  2. Hier ist eine Interpolation eine gerade Linie zwischen 2 Punkten. Wenn wir jedoch eine Koordinate eines Punktes auf dieser Linie kennen (y = 1,00), können wir einfach die andere Koordinate berechnen, indem wir die lineare Gleichung der Geraden auflösen, und das Problem wird in wenigen arithmetischen Operationen gelöst.

Unter Berücksichtigung dieser Gründe können wir eine effizientere Lösung für Ihr Problem entwickeln:

# solution of linear function between two points (2. reason)
def lin_itp(p1,p2,tv):
    """Get x coord of point on line

    Determine the x coord. of a point (x, target_value) on the line
    through the points p1, p2.

    Approach:
      - parametrize x, y between p1 and p2: 
          x = p1[0] + t*(p2[0]-p1[0])
          y = p1[1] + t*(p2[1]-p1[1])
      - set y = tv and resolve 2nd eqt for t
          t = (tv - p1[1]) / (p2[1] - p1[1])
      - replace t in 1st eqt with solution for t
          x = p1[0] + (tv - p1[1])*(p2[0] - p1[0])/(p2[1] - p1[1])
    """
    return float(p1[0] + (tv - p1[1])*(p2[0] - p1[0])/(p2[1] - p1[1])) 

# target value:
t_v = 1.0
for y in [0, 1]:
    arr_sd = arr.isel(y=y)
    # get index for the value closest to the target value (but smaller)
    s_udim = int(xr.where(arr_sd - t_v <=0, arr_sd, arr_sd.min()).argmax())
    # I'm explicitly defining the two points here
    ps_itp = arr_sd[s_udim:s_udim+2]
    p1, p2 = (ps_itp.x[0], ps_itp[0]), (ps_itp.x[1], ps_itp[1])
    print(lin_itp(p1,p2,t_v))

>>> 1.0344827586206897
>>> 1.3214285714285714

Jojo
quelle
1
Sie haben einen Fehler gemacht, als Sie sagten: "arr_sd = arr.isel (y = 0)" Sie meinen "arr_sd = arr.isel (y = y)"
Hoogendijk
@Hoogendijk du hast recht, danke. habe das nicht gesehen. Hoffe die Antwort war hilfreich. :)
Jojo
Ja, es war nützlich, aber ich entschied mich trotzdem zu prüfen, ob ich es verbessern und die Notwendigkeit einer for-Schleife beseitigen könnte.
Hoogendijk
0

Das Problem, das ich mit Jojos Antwort hatte, ist, dass es schwierig ist, es in vielen Dimensionen zu erweitern und die Röntgenstruktur beizubehalten. Daher habe ich mich entschlossen, dies weiter zu untersuchen. Ich habe einige Ideen aus Jojos Code verwendet, um die folgende Antwort zu geben.

Ich erstelle zwei Arrays, eines mit der Bedingung, dass die Werte kleiner sind als das, wonach ich suche, und eines mit der Bedingung, dass sie größer sein müssen. Ich verschiebe die zweite in x-Richtung um minus 1. Jetzt kombiniere ich sie in einer normalen linearen Interpolationsformel. Die beiden Arrays haben nur Werte, die sich am Rand der Bedingung überlappen. Wenn nicht um -1 verschoben, würden sich keine Werte überlappen. In der letzten Zeile summiere ich über die x-Richtung und da alle anderen Werte sind NaN, extrahiere ich den korrekten Wert und entferne dabei die x-Richtung aus dem DataArray.

def interpolate_dimension_x(arr, target_value, step):
    M0 = arr.where(arr - target_value <= 0)
    M1 = arr.where(arr - target_value > 0).shift(x=-1)

    work_mat = M0.x + step * (target_value - M0) / (M1 - M0)

    return work_mat.sum(dim='x')
interpolate_dimension_x(arr, 1, 0.25)

>>> <xarray.DataArray (y: 2)>
array([1.034483, 1.321429])
Coordinates:
  * y        (y) int32 1 2

Ich habe einige Nachteile mit meinem Code. Der Code funktioniert nur, wenn M0 und M1 einen Wert finden, der die Bedingung erfüllt. Andernfalls werden alle Werte in dieser Zeile auf gesetzt NaN. Um Probleme mit M0 zu vermeiden, habe ich beschlossen, die x-Werte nur bei 0 zu beginnen, da mein Zielwert immer größer als 0 ist. Um Probleme mit M1 zu vermeiden, wähle ich meine Werte für x groß genug, damit ich weiß, dass meine Werte dort sind . Dies sind natürlich keine idealen Lösungen und können den Code brechen. Wenn ich etwas mehr Erfahrung mit Xarray und Python habe, könnte ich umschreiben. Zusammenfassend habe ich folgende Punkte, die ich lösen möchte:

  • Wie extrapoliere ich Werte außerhalb des x-Bereichs? Ich stelle derzeit nur sicher, dass mein x-Bereich groß genug ist, damit die Antworten in diesen Bereich fallen.
  • Wie kann der Code für eine variable Schrittgröße robust gemacht werden?
  • Wie erstelle ich den Code, damit meine Dimension dynamisch ausgewählt werden kann (jetzt funktioniert es nur noch für 'x')?
  • Optimierungen sind willkommen.
Hoogendijk
quelle