`xarray-leaflet`#

xarray-leaflet ist eine xarray-Erweiterung für das Plotten von gekachelten Karten. Sowohl xarray als auch Leaflet können mit Datenfragmenten arbeiten, xarray durch Dask Chunks und Leaflet durch map tiles. Mit xarray-leaflet arbeiten beide zusammen.

Installation#

Ihr könnt xarray-leaflet in eurem Jupyter-Kernel installieren mit:

$ pipenv install xarray-leaflet
Installing xarray-leaflet…
…

Standardmäßig generiert xarray-leaflet Kacheln in temporären Verzeichnissen. Bei dynamischen Karten wird bei jeder Interaktion mit der Karte ein neues Verzeichnis erstellt, entweder durch Ziehen oder Zoomen. Dies liegt daran, dass eine direkte Zuordnung zwischen dem Kachelverzeichnis und der URL besteht, unter der die Kacheln bereitgestellt werden. Da bei dynamischen Karten Kacheln nicht vom Browser zwischengespeichert werden sollten, muss sich die URL ständig ändern. Diese temporären Verzeichnisse werden derzeit nicht automatisch bereinigt. Ihr solltet dies daher regelmäßig selbst tun. In ix-Systemen sind sie unter /tmp/xarray_leaflet_*.

Beispiel#

Um das Beispiel ausführen zu können, müsst Ihr zusätzlich die folgenden Pakete in Eurem Kernel installieren:

requests
tqdm
matplotlib
scipy

[1]:

import requests
import os
from tqdm import tqdm
import zipfile
import xarray as xr
import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt
import xarray_leaflet
from ipyleaflet import Map, basemaps

Nun laden wir den HydroSHEDS-Datensatz herunter:

[2]:

url = 'https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/hydrosheds/sa_30s_zip_grid/sa_acc_30s_grid.zip'
filename = os.path.basename(url)
name = filename[:filename.find('_grid')]
adffile = name + '/' + name + '/w001001.adf'

if not os.path.exists(adffile):
    r = requests.get(url, stream=True)
    with open(filename, 'wb') as f:
        total_length = int(r.headers.get('content-length'))
        for chunk in tqdm(r.iter_content(chunk_size=1024), total=(total_length/1024) + 1):
            if chunk:
                f.write(chunk)
                f.flush()
    zip = zipfile.ZipFile(filename)
    zip.extractall('.')

[3]:

da = xr.open_rasterio(adffile)
da

[3]:

xarray.DataArray

band: 1
y: 8520
x: 7320

...
```
[62366400 values with dtype=int32]
```

Coordinates: (3)

band
(band)
int64
1
```
array([1])
```

(y)

float64

15.0 14.99 14.98 ... -55.99 -56.0

array([ 14.995833,  14.9875  ,  14.979167, ..., -55.979167, -55.9875  ,
       -55.995833])

(x)

float64

-93.0 -92.99 ... -32.01 -32.0

array([-92.995833, -92.9875  , -92.979167, ..., -32.020833, -32.0125  ,
       -32.004167])

Attributes: (7)
transform :
(0.0083333333333333, 0.0, -93.0, 0.0, -0.0083333333333333, 14.999999999999716)
crs :
+init=epsg:4326
res :
(0.0083333333333333, 0.0083333333333333)
is_tiled :
0
nodatavals :
(-2147483647.0,)
scales :
(1.0,)
offsets :
(0.0,)

Ihr könnt sehen, dass die Projektion EPSG: 4326 (auch bekannt als WGS84) ist. Dies ist die einzige Projektion, die derzeit von xarray-leaflet unterstützt wird. Hier entspricht die Koordinate x dem Längengrad und y dem Breitengrad. Die band-Koordinate ist hier ziemlich nutzlos. Die einzige Vorverarbeitung, die wir durchführen werden, besteht darin, NaN-Werte durch Nullen zu ersetzen.

[4]:

nan = da.attrs['nodatavals'][0]
da = da.sel(band=1)
da = xr.where(da==nan, 0, da)

xarray-leaflet verfügt über zwei Modi, einen für statische Karten und einen für dynamische Karten. Mit dynamisch ist gemint, dass Ihr abhängig von der aktuellen Kartenansicht nicht dasselbe seht. Wenn Ihr beispielsweise zoomt, möchtet ihr vermutlich Details angezeigt bekommen, die vorher nicht zu sehen waren. Die Karte passt sich also wirklich an. Der Nachteil dynamischer Karten ist, dass Ihr bei der Interaktion mit der Karte mehr Flackern seht, da Kacheln aktualisiert werden müssen, sobald Ihr zieht oder zoomt.

Das folgende Beispiel könnt Ihr so konfigurieren, dass es sowohl im statischen als auch im dynamischen Modus funktioniert. Ihr könnt die folgende dynamische Variable von True auf False ändern, um in den statischen Modus zu wechseln.

[5]:

dynamic = True
if not dynamic:
    global_vmin = np.min(da).values
    global_vmax = np.max(da).values

Es gibt 3 Stufen in der array-leaflet, in denen ihr die Daten transofmieren könnt:

Der erste Stufe arbeitet am sichtbaren Teil. Wenn Ihr dynamische Transformationen anwenden wollt, könnt ihr diese hier als Python-Funktionen mit den Daten als erstem Argument und optional Parameter aus den vorherigen Transformationen angeben. Da transform0 die erste Transformation ist, werden DataArray keine weiteren Eingabeparameter zugewiesen. Hier geben wir die in den Daten vorhandenen Minimal- und Maximalwerte zurück:

[6]:

def transform0(da):
    if dynamic:
        vmin = np.min(da).values
        vmax = np.max(da).values
    else:
        vmin = global_vmin
        vmax = global_vmax
    return da, vmin, vmax

Die zweite Transformation gilt für die Daten, die in jeder Leaflet-Kachel vor der Neuprojektion enthalten sind. Bei der Neuprojektion müsst Ihr Daten in den Speicher laden. Daher möchtet Ihr die Daten vermutlich verkleinern auf ungefähr die gleiche Anzahl von Pixeln wie in einer Kachel (256 x 256). In dieser Transformation möchtet Ihr auch, dass Ihre Daten zwischen 0 und 1 normalisiert werden, was wir dank der vmin- und vmax-Werte tun können, die wir aus der vorherigen Transformation erhalten haben.

[7]:

def transform1(xr_tile, vmin, vmax):
    ny, nx = xr_tile.shape
    wx = nx // (256 // 2)
    wy = ny // (256 // 2)
    if wx > 0 and wy > 0:
        xr_tile = xr_tile.coarsen(x=wx, y=wy, boundary="trim").max()
    xr_tile = (xr_tile - vmin) / (vmax - vmin)
    return xr_tile

Bei der dritten und letzten Umwandlung werden die Daten in Web Mercator (die in Leaflet standardmäßig verwendete Projektion) transformiert. Dies ist die letzte Möglichkeit, Eure Daten zu transformieren, bevor Ihr sie in einer PNG-Datei speichert. Daher passen wir jetzt noch das Styling an und lassen die Flüsse dicker erscheinen, reduzieren die Amplitude mit einer Quadratwurzel und wählen eine Farbkarte.

[8]:

def transform2(np_tile):
    radius = 2
    circle = np.zeros((2*radius+1, 2*radius+1)).astype('uint8')
    y, x = np.ogrid[-radius:radius+1,-radius:radius+1]
    index = x**2 + y**2 <= radius**2
    circle[index] = 1
    np_tile = np.sqrt(np_tile)
    np_tile = scipy.ndimage.maximum_filter(np_tile, footprint=circle)
    np_tile = plt.cm.inferno(np_tile)
    return np_tile

Wir sind fast fertig, wir müssen nur eine Karte erstellen, bevor wir sie an unser DataArray übergeben:

[9]:

m = Map(center=[-20, -60], zoom=3, basemap=basemaps.CartoDB.DarkMatter)
m

Um unsere Daten auf der Karte anzuzeigen, rufen wir .leaflet.plot() in unserem DataArray auf und übergeben als Parameter die Karte, Längen- und Breitengrade, die Transformationsfunktionen und den dynamischen Wert. Wir erhalten eine Schicht zurück, die wir z – einen Schieberegler zum Einstellen der Deckkraft.

[11]:

l = da.leaflet.plot(m, lat_dim='y', lon_dim='x',
                    transform0=transform0,
                    transform1=transform1,
                    transform2=transform2,
                    dynamic=dynamic)
l.interact(opacity=(0.1,1.0,0.1))

xarray-leaflet#

Installation#

Beispiel#

`xarray-leaflet`#