Auswahl plotten (netcdf Dateien)

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Antworten
chulia
User
Beiträge: 29
Registriert: Mittwoch 17. August 2011, 10:26

Montag 28. April 2014, 09:04

Hallo ihr Lieben,
ich hab mal wieder ein Problem mit Python.
Ich arbeite gerade mit netcdf Dateien. Und zwar habe ich eine netcdf Datei in der für Deutschland die Bodenfeuchte über 60 Jahre zusammengefasst ist und in einer netcdf Datei habe ich die naturräumliche Gliederung Deutschlands (das ist quasi nur ein Rasterdatensatz, der auf die Bodenfeuchtedaten angepasst wurde).
Deutschland besteht aus 70 Hauptgruppen (durchnummeriert von 1-70).
Mein Ziel ist es, dass ich mir die Koordinaten einer Hauptgruppe herausfiltere und dieses Koordinaten nutze um mir die Bodenfeuchte für diese Hauptgruppe anzuzeigen. Leider kann ich nicht diese Auswahl plotten, es folgt immer die Fehlermeldung, dass z ein 2D Array sein muss.

Hier mein bisheriger Ansatz:

Code: Alles auswählen

# 1. netcdf Datei (Bodenfeuchte SMI):

f.dimensions # {'time': None, 'easting': 175, 'northing': 225}
# daraus herausgefiltert, SMI mit Zeit, Northing und Easting
SMI = f.variables['SMI']
SMI.dimensions #('time', 'northing', 'easting')
SMI.shape #(720, 225, 175)

#2. netcdf Datei (Naturraum):
Natraum.dimensions #{'y': 225, 'x': 175}

Latitude=Natraum.variables['x'][:]
Longitude=Natraum.variables['y'][:]
Hauptgruppe=Natraum.variables['Haupteinheit'][:] 
Hauptgruppe.shape # (225, 175)

# ob Longitude und Nothing und Latitude und Easting der beiden netcdf gleich sind habe ich über eine Schleife geprüft

# jetzt filtere ich die Hauptgruppe 1 heraus:

b=np.where(Hauptgruppe==1)
c=b[0] 
d=b[1]

# Gegenkontrolle:

Hauptgruppe[c[:],d[:]] # ergibt ein Array voller 1en

# so und nun würde ich das gerne plotten (sodass ich auch weiß, wo die Hauptgruppe liegt)

# alle Hauptgruppen konnte ich mir darüber anzeigen lassen:

levels=MaxNLocator(nbins=70).bin_boundaries(1,70)
plt.contour(Latitude,Longitude,Hauptgruppe,levels=levels)

# und nun dachte ich geb ich einfach dafür die Indexe von c und d ein.

plt.contour(Latitude[c[:]],Longitude[d[:]],Hauptgruppe[c[:],d[:]])

# aber da kommt folgende Fehlermeldung:
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/.../lib/python2.7/site-packages/matplotlib/pyplot.py", line 2632, in contour
    ret = ax.contour(*args, **kwargs)
  File "/.../lib/python2.7/site-packages/matplotlib/axes.py", line 7976, in contour
    return mcontour.QuadContourSet(self, *args, **kwargs)
  File "/.../lib/python2.7/site-packages/matplotlib/contour.py", line 1414, in __init__
    ContourSet.__init__(self, ax, *args, **kwargs)
  File "/.../lib/python2.7/site-packages/matplotlib/contour.py", line 860, in __init__
    self._process_args(*args, **kwargs)
  File "/.../lib/python2.7/site-packages/matplotlib/contour.py", line 1427, in _process_args
    x, y, z = self._contour_args(args, kwargs)
  File "/.../lib/python2.7/site-packages/matplotlib/contour.py", line 1491, in _contour_args
    x, y, z = self._check_xyz(args[:3], kwargs)
  File "/.../lib/python2.7/site-packages/matplotlib/contour.py", line 1525, in _check_xyz
    raise TypeError("Input z must be a 2D array.")
TypeError: Input z must be a 2D array.

# eigentlich will ich mir ja nicht nur die Hauptgruppe anzeigen lassen, sondern den SMI, aber auch da funktioniert es nicht

plt.contour(easting[c[:]],northing[d[:]],SMI[:,c[:],d[:]],MaxNLocator(nbins=15).bin_boundaries(0,1))
# als Fehlermeldung kommt dann:
TypeError: Length of y must be number of rows in z.

Ich kann mir auf das ganze keinen Reim machen, was mache ich falsch?
Kann mir einer helfen?

Vielen Dank =)
BlackJack

Montag 28. April 2014, 10:42

@chulia: Die Fehlermeldung ist doch recht deutlich. Das dritte Argument muss ein 2D-Array sein. Und Du fütterst da ein eindimensionales Array rein wobei jeder Wert davon auch noch 1 ist. Wie hattest Du Dir denn vorgestellt das man *daraus* einen Kontourplot erstellen kann?

Man könnte zum Beispiel versuchen das Ergebnis von ``Hauptgruppe == 1`` zu verwenden, was ein 2D-Array ist wo ein `True` steht wenn der Wert in `Hauptgruppe` an der Stelle 1 ist und `False` wenn dort eine andere Zahl steht.

Du erstellst übrigens ziemlich inflationär und unnötig dauernd neue Views in die Arrays. Warum?
chulia
User
Beiträge: 29
Registriert: Mittwoch 17. August 2011, 10:26

Montag 28. April 2014, 10:56

Danke für deine Hilfe.

Ich bin blutige Anfängerin in Python und versuch mir das gerade selber beizubringen.
Ich will mir ja einfach nur ansehen, wo in Deutschland die Hauptgruppe 1 ist. Deswegen der Contourplot. Im Grunde will ich damit nur erreichen, dass es funktioniert und ich das dann auf mein SMI (dort stehen dann ja unterschiedliche werte drin (nicht nur 1)) anwenden kann.
Wenn ich das mit meinem b versuche also so:

Code: Alles auswählen

plt.contour(Latitude[c[:]],Longitude[d[:]],b[:]),
dann kommt die Fehlermeldung, dass Length of y must be number of rows in z... wobei y 279 lang ist und b ein Tuple aus 2x 279 (was ja eigentlich stimmen sollte ...)

und was meinst du mit
Du erstellst übrigens ziemlich inflationär und unnötig dauernd neue Views in die Arrays. Warum?
Ich fänds super, wenn einer etwas Licht ins Dunkel bringen kann ^^
BlackJack

Montag 28. April 2014, 11:47

@chulia: Wenn das stimmen würde, gäbe es die Fehlermeldung nicht. Und Du sagst ja selber das die Anzahl der Zeilen in `b` *zwei* ist. Und die Koordinaten machen als `z` keinen Sinn, denn wie schon mal gesagt wird dort ein 2D-Array mit Werten erwartet aufgrund derer ein Kontourplot erstellt wird.

Bezüglich der unnötigen Views: Was sollen die ganzen ``[:]``? Das macht keinen Sinn.
chulia
User
Beiträge: 29
Registriert: Mittwoch 17. August 2011, 10:26

Montag 28. April 2014, 12:28

mmh zu den [:], dachte ich, dass man das immer braucht. (Ich komm aus der Matlabwelt ^^, aber hab eben gesehen, dass es auch ohne geht)

Okay, dass das mit meinen Hauptgruppen nicht geht, hab ich ja verstanden, was ich nicht allerdings nicht verstehe ist, warum bei SMI das gleiche Problem vorliegt.

wenn ich den normalen SMI mir ansehe, dann ist das ein 3D array

Code: Alles auswählen

SMI.shape
(720, 225, 175)
720 Zeitschritte, 225 xKoordinaten und 175 yKoordinaten und wenn ich für x und y etwas eingebe, erhalte ich den dazugehörigen SMI Wert:

Code: Alles auswählen

SMI[5,44,123] # 5. Zeitschritt, Koordinatenpaar (44,123)
0.88106465 # Bodenfeuchte
Und wenn ich für Latitude und Longitude jeweils 44 und 123 eingebe, komm ich auch auf richtige Koordinaten

Code: Alles auswählen

Latitude[44]
4216000.0
Longitude[123]
5626000.0
deswegen müssen beim plot Befehl immer die Koordinaten mit rein:

Code: Alles auswählen

plt.contour(Latitude[c],Longitude[d],SMI[5,c,d],MaxNLocator(nbins=15).bin_boundaries(0,1))
Lass mich aber jetzt die Dimensionen von SMI anzeigen, ist aus dem 3D Array ein 1D Array geworden und deswegen kommt wieder die Fehlermeldung, dass z ein 2D array sein muss.

Code: Alles auswählen

SMI[5,c,d].shape
(279,)
Wie geht das? Also irgendwie ist mir schon klar, wie das geht, weil ich ja für den 5. Zeitschritt mir die SMI Werte für die Koordinatenpaare c und d ausgeben lasse.
Also ist es nunmal 1D und nicht 2D.

Gibt es vielleicht einen anderen plotbefehl mitdem ich das einfacherer darstellen kann?

Also gebe ich ein

Code: Alles auswählen

plt.contour(SMI[5,:,:])
gibt es kein Problem, aber sobald ich das 2. und 3. Array spezifiziere und eben nur eine Auswahl geplottet werden soll, kommt die 2D Array Fehlermeldung


Danke dir =)
BlackJack

Montag 28. April 2014, 14:28

@chulia: Ob ein Kontourplot das richtige für Deine Daten ist hängt letztendlich von der Bedeutung der Daten ab und was Du in so einem Plot zeigen möchtest.

Du bereitest die Daten falsch auf. Der Ansatz mit den Koordinaten ist Unsinn, denn Du brauchst letztendlich die Daten in der Form wie sie in einem Zeitschnitt aus dem SMI-Array kommen.

Du könntest Dir ja erst einmal einen Kontourplot von den ganzen Daten ansehen.

Wenn Du dann einen nur von einem Bereich erstellen möchtest, stellen sich ein paar Fragen. Man könnte über das Hauptbereiche-Array herausfinden welche Werte zu einem Bereich gehören anhand der Daten dann aus den 2D-Daten zum Plotten einen 2D-Unterauschnitt bilden. Und vielleicht die Daten die nicht zum Hauptbereich in diesen Daten auf NaN setzen.
chulia
User
Beiträge: 29
Registriert: Mittwoch 17. August 2011, 10:26

Montag 28. April 2014, 14:35

Ja genau, das scheint ja mein Problem zu sein. Wie bekomme ich denn einen 2D Unterausschnitt hin?
Deswegen bin ich ja über die Koordinaten gegangen.
Kannst du vielleicht mal ein Code Beispiel schicken? Ich glaub dann kann ich mir das etwas besser vorstellen.

Danke =)
BlackJack

Montag 28. April 2014, 15:31

@chulia: Also aus einem 2D-Array ein 2D-Array ”ausschneiden” geht ganz normal über „slicing”. Das ist eine Grundoperation die einem eigentlich ein Grundlagentutorial zu Numpy beinbringen sollte. Es gibt in der Numpy-Dokumentation ein Tutorial.

Letztendlich ist es „slicing” um die Zeilen und Spalten in der jeweiligen Dimension auszuwählen:

Code: Alles auswählen

In [90]: A
Out[90]: 
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])

In [91]: A[2:4]
Out[91]: 
array([[10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19]])

In [92]: A[2:4, 1:3]
Out[92]: 
array([[11, 12],
       [16, 17]])
Bei [91] werden die Zeilen 3 bis 4 ausgewählt (dran denken das Indizes in Python bei 0 anfangen und der Endpunkt bei Slices nicht enthalten ist), und bei [92] werden in der zweiten Dimension noch die Spalten ausgewählt.

Bleibt noch die Frage wie man die Eckkoordinaten für die Bereiche ermittelt. Da braucht man tatsächlich die Koordinaten, also `numpy.argwhere()` auf eine Bitmap angewendet, die nur einen Bereich enthält, und davon dann das Minimum und das Maximum über die erste Achse. Ungetestet:

Code: Alles auswählen

coordinates = numpy.argwhere(Hauptgruppe == 1)
(x_min, y_min), (x_max, y_max) = coordinates.min(0), coordinates.max(0) + 1
Damit hätte man die passenden Werte die man zum „slicen” der Daten für die Hauptgruppe 1 verwenden kann.

An der Stelle sollte man sich vielleicht auch mal Gedanken um Funktionen machen, denn das alles so auf Modulebene runter zu hacken scheint bei Naturwissenschaftlern zwar beliebt zu sein, aber spätestens nach einem Jahr steigen da selbst die Autoren oft nicht mehr durch und müssen sich wieder mühsam in ihren eigenen Quelltext einarbeiten.
chulia
User
Beiträge: 29
Registriert: Mittwoch 17. August 2011, 10:26

Montag 28. April 2014, 15:43

BlackJack, ich danke dir erstmal für deine schnelle Hilfe.
Ich versuchs morgen wieder, wenn ich auf Arbeit bin.
Zumindestens kann ich mir jetzt etwas drunter vorstellen.
Vielen Dank =)
Antworten