Karten und Höhendiagramme

Neben Häufigkeit und Verteilung unterscheidet sich auch die Datengrundlage von Art zu Art. Von einigen Arten kennen wir sämtliche Neststandorte, bei anderen können wir basierend auf den Kartierungen die Dichte pro Kilometerquadrat berechnen. Bei schwer erfassbaren Arten können wir lediglich Aussagen zur Vorkommenswahrscheinlichkeit machen, kennen die Dichte aber nicht. Um diesen Unterschieden Rechnung zu tragen, entwarfen wir verschiedene Kartentypen für die räumliche Darstellung des Vorkommens und dessen Veränderung seit 1993–1996. Wir entschieden für jede Art, welcher Kartentyp am aussagekräftigsten ist. In der Tabelle im Abschnitt «Erfassungsmethoden der Arten» (S. 53–59) ist angegeben, mit welchem Kartentyp der Zustand 2013–2016 und die Veränderung des Vorkommens seit 1993–1996 abgebildet werden.

Dichtekarte 2013–2016

Für 75 relativ weit verbreitete und zumindest lokal häufige Arten konnten wir Dichtekarten erstellen. Dafür verwendeten wir die Kartierungsdaten.

Basierend auf der Anzahl gefundener Reviere pro Kartierungsrundgang und mit 16 verschiedenen Umweltvariablen schätzten wir die Revierzahl pro Art für alle Kilometerquadrate der Schweiz und Liechtensteins. Dazu verwendeten wir ein Binomial Mixture Model. Dies entspricht einer Poisson-Regression, bei der zusätzlich die Entdeckungswahrscheinlichkeit der Vogelart miteinbezogen wird. So berücksichtigten wir, dass bei einem Kartierungsrundgang jeweils nur ein Teil der tatsächlich anwesenden Vögel entdeckt wurde.

Wir ermittelten die Werte der 16 Umweltvariablen für sämtliche ganz oder teilweise in der Schweiz gelegene Kilometerquadrate sowie für ausländische Kilometerquadrate, in denen Kartierungen stattfanden. Für sieben flächige Variablen (Nr. 1–7 in der Tabelle Umweltvariablen) wurde der flächenmässige Anteil für jedes Kilometerquadrat der Schweiz und Liechtensteins aus der Arealstatistik nach der Nomenklatur NOAS04 eruiert. Für die Jahre 2013–2016 verwendeten wir die Arealstatistik 2004–2009. Für Kilometerquadrate im Ausland (ohne Liechtenstein) wurden die Werte möglichst analog zur Arealstatistik aus den CORINE Land Cover Daten 2012 (Datenaufnahme 2011–2012) abgeleitet. Fünf weitere Variablen (Nr. 8–12 in der Tabelle Umweltvariablen) berechneten wir für die Periode 2013–2016 aus dem Topografischen Landschaftsmodell oder dem VECTOR25 Datensatz. Für wenige Quadratkilometer im Ausland wurden die entsprechenden Grundlagendaten auf Basis der aktuellen Landeskarten 1 : 25 000 nachdigitalisiert. Die Angaben zum Stickstoffeintrag (Nr. 13 in der Tabelle Umweltvariablen) ermittelte die Firma Meteotest im Auftrag des Bundesamtes für Umwelt; sie gelten für das Jahr 2010. Drei Variablen (Nr. 14–16 in der Tabelle Umweltvariablen) basieren auf dem digitalen Höhenmodell. Das digitale Höhenmodell überzieht die Schweiz mit einem Netz mit Maschenweite 25 m. Für jeden Knoten ist ein Wert (z.B. Meereshöhe) angegeben. Für jedes Kilometerquadrat mittelten wir die Werte aller innerhalb des Quadrats liegender Knotenpunkte und verwendeten diesen Mittelwert als Kovariablenwert.

Oft war der Einfluss auf die Dichte einer Art nicht für alle Variablen schätzbar, insbesondere da gewisse Variablen recht stark korreliert waren. Daher wurden Variablen für jede Art schrittweise weggelassen, wenn ihr Einfluss nur schwer schätzbar war und gemäss der Biologie einer Art kein Einfluss auf die Dichte erwartet wurde. Wir verwendeten Penalized 2D-Splines, um räumliche Autokorrelation explizit modellieren zu können. Dies war wichtig, da eine Art trotz passendem Lebensraum regionale Unterschiede in der Dichte zeigen kann. Durch den Miteinbezug von grenznahen Kilometerquadraten im Ausland liessen sich Randeffekte vermeiden, die bei Auswertungen mit Splines auftreten können. Schliesslich setzten wir die Dichte in sämtlichen Kilometerquadraten oberhalb einer für jede Art festgelegten Meereshöhe auf null, um Artefakte im Hochgebirge zu vermeiden. Mit den vom Modell geschätzten Parameterwerten konnten wir so für jedes der rund 41 000 Kilometerquadrate der Schweiz und Liechtenstein eine Schätzung der Revierzahl machen. Bei den meisten Arten waren die Schätzungen in einigen wenigen Kilometerquadraten unrealistisch hoch. Für die optische Darstellung korrigierten wir diese Werte nach unten. Alle Werte, die über einer festgelegten Obergrenze lagen, wurden auf den Wert dieser Obergrenze reduziert. Für fünf Arten setzten wir die Obergrenze von Hand fest (Waldlaubsänger, Fitis, Haus- und Italiensperling, Gebirgsstelze). Für die übrigen Arten wählten wir als Obergrenzen das 99,5-%-Quantil aus den Schätzwerten aller Kilometerquadrate. Quadrate mit einer geschätzten Dichte von weniger als 0,05 wurden für die Ermittlung des 99,5-%-Quantils vorgängig weggelassen. Für eine optisch ansprechende Darstellung wurden die Schätzwerte mit einer Interpolation leicht geglättet.

Dichteänderungskarte

Für 70 der 75 Arten, für die wir eine Dichtekarte für 2013–2016 berechneten, konnten wir auch eine Dichteänderungskarte erstellen, welche die Veränderung in der Dichte zwischen 1993–1996 und 2013–2016 aufzeigt. Im Gegensatz zu den Daten 2013–2016 war es mit den Daten 1993–1996 nicht möglich, die Entdeckungswahrscheinlichkeit bei der Schätzung der Dichten zu berücksichtigen. Dies, weil die Daten 1993–1996 nur die Anzahl der gefundenen Reviere pro Art und Kilometerquadrat enthielten, nicht aber die Anzahl der Reviere, die pro Kartierungsrundgang durch eine Beobachtung belegt waren («Reviere pro Rundgang»). Analog zur Erstellung der Dichtekarten schätzten wir für jedes Kilometerquadrat eine Dichte für die Periode 1993–1996 sowie für die Periode 2013–2016, diesmal jedoch ohne Berücksichtigung der Entdeckungswahrscheinlichkeit. Als Datengrundlage dienten die Kartierungen aus den jeweiligen Perioden. Im Gegensatz zum Modell mit Berücksichtigung der Entdeckungswahrscheinlichkeit, in dem als Zielvariable die Anzahl «Reviere pro Rundgang» modelliert wird, verwendeten wir hier als Zielvariable die Anzahl gefundener Reviere pro Kilometerquadrat.

Als erklärende Variablen verwendeten wir analog zu den Dichtekarten dieselben 16 Umweltvariablen. Bei den Variablen Flüsse und Bäche, Seeufer, Meereshöhe, Exposition und Hangneigung verwendeten wir für beide Perioden dieselben Variablenwerte. Bei den übrigen Variablen ermittelten wir die Werte für 1993–1996 separat, basierend auf der Arealstatistik 1992–1997, den CORINE Land Cover Daten 1990 sowie dem VECTOR25-Datensatz. Die Angaben zum Stickstoffeintrag erhielten wir von der Firma Meteotest; sie gelten für das Jahr 1990. Pro Art benützten wir für 1993–1996 und 2013–2016 dasselbe Set an Umweltvariablen. Für die Auswertung verwendeten wir ein Generalized Linear Model anstatt eines Binomial Mixture Models.

Weiter mussten wir beachten, dass 1993–1996 für alle Arten eine Obergrenze bei der Anzahl Reviere festgelegt worden war. Hatte ein Bearbeiter beispielsweise mehr als zehn revieranzeigende Buchfinken gefunden, musste er diese Art nicht weiter kartieren. Um trotzdem die tatsächliche Revierzahl schätzen zu können, wurde diese sogenannte Zensurierung beim Formulieren des Modells miteinbezogen. Um die Vergleichbarkeit der Daten zu erreichen, zensurierten wir die Daten der 2013–2016 durchgeführten Kartierungen analog zum Vorgehen 1993–1996. Teilweise waren die so ermittelten Ergebnisse unbefriedigend und stimmten nicht mit der Einschätzung von Experten oder mit regionalen, aus dem Projekt «Monitoring Häufige Brutvögel» (MHB) geschätzten Trends überein. In diesen Fällen ergänzten wir die zensurierten Daten der Kartierungen zusätzlich mit unzensurierten Daten aus dem MHB. Die Kartierungsdaten der Periode 1993–1996 ergänzten wir mit den MHB-Daten aus dem Jahr 2000, die Daten der Periode 2013–2016 mit den MHB-Daten 2016. So konnten wir in einigen Fällen die Schätzung der Dichte und der Dichteänderung verbessern.

Für die Darstellung glätteten wir die Dichteschätzungen für 1993–1996 und 2013–2016 separat. In einem ersten Schritt berechneten wir für jedes Kilometerquadrat den Durchschnitt der Dichteschätzungen aller neun Quadrate, die in einer 3 × 3-km-Matrix mit dem entsprechenden Quadrat als Mittelpunkt enthalten sind und setzten diesen Wert für das zentrale Kilometerquadrat ein. Schliesslich berechneten wir für jedes Kilometerquadrat die Differenz dieser geglätteten Werte zwischen den zwei Perioden 1993–1996 und 2013–2016.

Verbreitungsänderungskarte

Für 21 meist eher seltene Arten konnten wir Verbreitungsänderungskarten erstellen. Dazu schätzten wir die artspezifische Vorkommenswahrscheinlichkeit pro Kilometerquadrat zusätzlich auch für die Periode 1993–1996. Wir verwendeten dafür die Daten der Kartierungen, der vollständigen Beobachtungslisten sowie die Einzelmeldungen der Jahre 1993–1996 und werteten diese mit einem Site-Occupancy Model aus.

Bei 608 der 2934 Kartierungen von 1993–1996 wurden die auf Papier vorliegenden Kartierungsergebnisse für alle Rundgänge nachdigitalisiert. Für diese Flächen waren somit wiederholte Meldungen innerhalb derselben Saison verfügbar, was uns erlaubte, mittels Site-Occupancy Model die Entdeckungswahrscheinlichkeit zu schätzen. Für einen grossen Teil der Kartierungen von 1993–1996 hingegen war nur bekannt, ob eine Art irgendwann während den drei Rundgängen nachgewiesen wurde oder nicht. Um diesen Datentyp ebenfalls auswerten zu können, erweiterten wir das klassische Site-Occupancy Model folgendermassen: Wenn eine Art gar nicht gefunden wurde, war sie entweder wirklich nicht vorhanden, oder sie wurde drei Mal übersehen. Dies kommt mit einer Wahrscheinlichkeit von (1 – p) vor, wobei p die Entdeckungswahrscheinlichkeit der Art während eines Rundgangs ist, die basierend auf den Daten der 608 nacherfassten Kartierungen geschätzt werden konnte. So konnten wir basierend auf den verschiedenen vorhandenen Datentypen auch für die Jahre 1993–1996 eine Vorkommenswahrscheinlichkeit pro Kilometerquadrat ermitteln.

Analog zum Vorgehen bei den Dichteänderungskarten glätteten wir die Schätzungen der Vorkommenswahrscheinlichkeit mittels 3 × 3-km-Matrix und berechneten die Differenz dieser geglätteten Werte zwischen den zwei Perioden 1993–1996 und 2013–2016.

Verbreitungskarte basierend auf Dichteschätzung

Für verbreitete Arten mit grossen Revieren lagen oft genügend Meldungen aus den Kartierungen der Kilometerquadrate vor, um eine Dichtekarte zu erstellen. Allerdings resultierte so meistens eine zu hohe Dichte, da beim Kartieren oft auch Individuen mitgezählt wurden, deren Reviere mehrheitlich ausserhalb des Kilometerquadrats lagen. Für Birkhuhn, Schwarzmilan, Mäusebussard, Schwarzspecht, Turmfalke und Kolkrabe entschieden wir uns daher, die Dichteschätzungen mittels der Formel v = 1 – e von Whright in eine Schätzung der Vorkommenswahrscheinlichkeit zu reduzieren. Dabei entspricht N der geschätzten Revierzahl pro Kilometerquadrat und v der daraus abgeleiteten Vorkommenswahrscheinlichkeit. Für diese Arten bildeten wir den Zustand 2013–2016 und die Veränderung des Vorkommens seit 1993–1996 basierend auf den so transformierten Daten ab.

Verbreitungskarte 2013–2016

Für 47 Arten erstellten wir Verbreitungskarten, auf denen wir anstatt der Dichte die Vorkommenswahrscheinlichkeit pro Kilometerquadrat abbildeten. Es handelt sich hierbei um eher seltene Arten, für die aufgrund der zu geringen Zahl an Nachweisen in den Kartierungen der Kilometerquadrate keine Dichtekarten erstellt werden konnten. Als Datenquellen verwendeten wir dabei einerseits die Daten aus den Kartierungen, andererseits auch die vollständigen Beobachtungslisten und die Einzelmeldungen der bei der Vogelwarte als freiwillige Mitarbeitende eingeschriebenen Beobachterinnen und Beobachtern (Ornithologischer Informationsdienst ID). Da wir nicht die Dichte, sondern nur das Vorkommen (Präsenz/Absenz) modellierten, reduzierten wir die Kartierungsdaten zu Präsenz/Absenz-Daten (Werte >1 wurden gleich 1 gesetzt). Wir berücksichtigten also nur, ob eine Art in einem Rundgang gefunden wurde oder nicht. Auch aus den vollständigen Beobachtungslisten erhielten wir für jede Art direkt eine Angabe zu Präsenz/Absenz. Für die Einzelmeldungen konstruierten wir Absenzmeldungen analog Kéry et al.: Meldete ein ID-Mitarbeiter beispielsweise einen Sperber, aber keinen Mittelspecht, so werteten wir dies als Absenzmeldung für den Mittelspecht. Nur die Arten der Meldekategorien A und B konnten so behandelt werden, bei Arten der Meldekategorie C mussten wir uns auf die Daten aus den Kartierungen und den vollständigen Beobachtungslisten beschränken (für die Definition der Meldekategorien s. Schmid et al.). Durch die Verwendung der drei Datentypen (Kartierungen, vollständige Beobachtungslisten und Einzelmeldungen) konnten die Vorteile aus den verschiedenen Datensätzen kombiniert werden: Die Kartierungen sorgten für eine gute geografische Verteilung und deckten die verschiedenen Lebensräume ab, die Beobachtungslisten und die Einzelmeldungen erhöhten die räumliche Abdeckung um ein Vielfaches.

Für sieben Feuchtgebietsarten (Zwergtaucher, Wasserralle, Teichhuhn, Zwergdommel, Drosselrohrsänger, Rohr- und Feldschwirl) berücksichtigten wir zudem die Daten aus den einzelnen Kartierungsrundgängen des Projekts «Monitoring Brutvögel in Feuchtgebieten» (MF) als Angaben zu Präsenz/Absenz. Dies, weil einige besonders wichtige Feuchtgebiete dank ihres Schutzstatus im Frühling kaum begangen werden und daher mit den übrigen Datenquellen schlecht abgedeckt sind.

Als erklärende Variablen verwendeten wir dieselben 16 Umweltvariablen wie für die Dichtekarten. Analog zum Vorgehen bei den Dichtekarten reduzierten wir das Umweltvariablenset für gewisse Arten, modellierten die räumliche Autokorrelation mit Penalized 2D-Splines und setzten das Vorkommen in sämtlichen Kilometerquadraten oberhalb einer artspezifisch festgelegten Meereshöhe auf null. Als Analysemethode wählten wir ein Site-Occupancy Model. Analog zum Binomial Mixture Model konnte so die Entdeckungswahrscheinlichkeit miteinbezogen werden. Für jeden Datentyp (Kartierungen, vollständige Beobachtungslisten und Einzelmeldungen) wurde die Entdeckungswahrscheinlichkeit separat geschätzt. Analog zum Vorgehen bei den Dichtekarten verwendeten wir als räumliche Einheit für die Modellierung ein Kilometerquadrat. Wir schätzten also für alle Kilometerquadrate der Schweiz, mit welcher Wahrscheinlichkeit eine Art dort tatsächlich als Brutvogel vorkam.

Abgebildet haben wir in den meisten Fällen die vom Modell geschätzte Vorkommenswahrscheinlichkeit. In einigen Fällen bildeten wir das realisierte Vorkommen ab. Das realisierte Vorkommen ist eine Mischung von Modellschätzung und Rohdaten. Konkret bedeutet dies, dass auf einer Fläche mit einer Meldung der Art immer ein sicheres Vorkommen abgebildet wurde, unabhängig davon, ob das Modell für die Fläche eine hohe Vorkommenswahrscheinlichkeit schätzte oder nicht. Wurden auf einer Fläche hingegen keine Meldungen gemacht, dann entsprach das realisierte Vorkommen der vom Modell geschätzten Vorkommenswahrscheinlichkeit. Einerseits wurden Arten wie die Heidelerche so behandelt, bei denen ein grosser Teil der Vorkommen von den Beobachtern auch tatsächlich gefunden und gemeldet wurde. Andererseits gingen wir mit verschiedenen Feuchtgebietsarten so vor, weil viele Feuchtgebiete von Ornithologen vergleichsweise stark frequentiert werden oder die Bearbeitung via «Monitoring Brutvögel in Feuchtgebieten» (MF) recht vollständig ist und die dort vorkommenden Arten daher mit grosser Wahrscheinlichkeit gefunden werden.

Für eine optisch ansprechende Darstellung wurden die Werte analog dem für die Dichtekarten gewählten Vorgehen leicht geglättet.

Punktkarte 2013–2016

Die Verbreitung von 75 seltenen Arten und Koloniebrütern illustrierten wir mit einer Punktkarte. Als Datengrundlage verwendeten wir die Meldungen aus den Kartierungen der Kilometerquadrate, vollständige Beobachtungslisten und Einzelmeldungen sowie Spezialerhebungen (z.B. Kartierungen für das «Monitoring Brutvögel in Feuchtgebieten» (MF), Saatkrähen-Monitoring). Daraus ermittelten wir pro Ort und Jahr die Anzahl Reviere bzw. Brutpaare. Auf der Karte entspricht ein Punkt jeweils der Anzahl Reviere bzw. Brutpaare pro räumliche Einheit, gemittelt über die Jahre 2013–2016. Als räumliche Einheit wählten wir bei Koloniebrütern eine Kolonie. In Siedlungen haben wir die Daten üblicherweise für ganze Städte oder Stadtteile zusammengefasst. Bei Wasservögeln waren die räumlichen Einheiten einzelne Gewässer oder Gewässerabschnitte. Für sehr seltene Arten setzten wir einen Punkt, wenn an einem Ort zumindest in einem der Jahre 2013–2016 ein Nachweis mit ausreichend hohem Atlascode gemacht werden konnte.

Punktveränderungskarte

Für 21 Arten erstellten wir Veränderungskarten auf Punktniveau. Dafür bildeten wir die Veränderung der Anzahl Reviere bzw. Brutpaare pro räumliche Einheit ab. Dies war nur für Arten möglich, bei denen Ort und mittlere Anzahl Brutpaare der einzelnen Kolonien bzw. Brutstandorte bereits in den Jahren 1993–1996 bekannt waren. Oft war der Ort der Vorkommen 1993–1996 mit einer gewissen Unsicherheit behaftet. Nahe gelegene Standorte wurden auf der Punktveränderungskarte zusammengefasst, um dieser Unsicherheit Rechnung zu tragen und um die Karten übersichtlich zu halten. Gewisse auf der Punktkarte separat abgebildete Orte sind daher auf der Punktveränderungskarte nicht einzeln aufgeführt.

Höhendiagramm

Nebst der räumlichen Verteilung wird auch die Verteilung der einzelnen Arten über den Höhengradienten gezeigt. Dazu haben wir die auf den Karten abgebildeten Daten in Abschnitte von 100 Höhenmetern unterteilt. Wir summierten das geschätzte Vorkommen pro Höhenstufe und teilten es dann durch die für die gesamte Schweiz ermittelte Summe des geschätzten Vorkommens für die betreffende Art. So konnten wir für jede Höhenstufe den Anteil an der Gesamtpopulation ermitteln. Diesen Anteil bildeten wir im Höhendiagramm ab.

Für viele Arten konnten wir zusätzlich die Veränderung der Höhenverbreitung darstellen. Dazu berechneten wir für jede Höhenstufe, wie sich der Bestand zwischen 1993–1996 und 2013–2016 verändert hat. Diesen Wert teilten wir dann durch die Summe des aktuellen Bestands (2013–2016). Daraus ergab sich eine relative Veränderung des Bestands pro Höhenstufe als Anteil des Gesamtbestands der Jahre 2013–2016. Diesen Wert bildeten wir auf der rechten Seite des Höhendiagramms ab. Durch dieses Vorgehen ist die Länge der Balken zwischen den zwei Grafiken direkt vergleichbar. Die Höhenverteilung für die Periode 1993–1996 kann einfach ausgelesen werden, indem die Zunahmen vom aktuellen Zustand abgezogen und die Abnahmen dazugezählt werden.

Die Veränderung der Höhenverbreitung konnte nur für Arten dargestellt werden, bei denen die Datengrundlage bereits im Atlas 1993–1996 sehr gut war. In einigen Fällen mussten wir uns daher darauf beschränken, lediglich den aktuellen Zustand der Höhenverbreitung auszuweisen.

Validierung

Die Karten und Höhendiagramme wurden zur Validierung verschiedenen Artexperten vorgelegt. War das Ergebnis nicht zufriedenstellend, versuchten wir die Schätzungen durch eine sinnvolle Auswahl der Umweltvariablen zu verbessern. Wir orientierten uns dabei besonders an den ökologischen Ansprüchen der Art. Liess sich auf diese Weise keine Optimierung erzielen, wurde ein konservativeres Vorgehen gewählt. Konnte für eine Art keine befriedigende Dichtekarte erstellt werden, erstellten wir eine Verbreitungskarte. War die Verbreitungskarte ungenügend, verwendeten wir die Rohdaten und stellten die Verbreitung auf Stufe Atlasquadrat dar. Auf diese Weise gingen wir sowohl für die Abbildung des Zustands 2013–2016 als auch für die Abbildungen der Veränderung seit 1993–1996 vor.

Oft war es für die Modelle schwierig, die Obergrenze der Verbreitung einer Art korrekt zu schätzen, da die Landesfläche und damit auch die Menge der vorhandenen Daten gegen oben stark abnimmt. Nicht zutreffende Schätzungen werden insbesondere im Höhendiagramm sichtbar. In solchen Fällen passten wir die Höhenobergrenze an, oberhalb der die Dichte bzw. die Vorkommenswahrscheinlichkeit manuell auf null gesetzt wurden.

Verwendete Software

Wir bereiteten die Daten in R 3.3.2 auf. Für die statistische Auswertung verwendeten wir JAGS 4.2.x. Die Daten für die Punktkarten wurden in QGIS zusammengestellt.