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]:
- band: 1
- y: 8520
- x: 7320
- ...
[62366400 values with dtype=int32]
- band(band)int641
array([1])
- y(y)float6415.0 14.99 14.98 ... -55.99 -56.0
array([ 14.995833, 14.9875 , 14.979167, ..., -55.979167, -55.9875 , -55.995833])
- x(x)float64-93.0 -92.99 ... -32.01 -32.0
array([-92.995833, -92.9875 , -92.979167, ..., -32.020833, -32.0125 , -32.004167])
- 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, werdenDataArray
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
und1
normalisiert werden, was wir dank dervmin
- undvmax
-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))