Dynamically Modify Xarray Values

[jsimkins2]

Hello,

I’d like to load in a spatial grid as an xarray object, plot the grid, and modify its values interactively. I know this can be done with xarray using select features, but can it be done interactively? Specifically, can this be done with hvPlot or another package of the holoviz family?

import xarray as xr #version 0.15.1
import hvplot.xarray #version 0.6.0
import geoviews as gv #version 1.8.1
gv.extension('bokeh')

# read in downloaded dataset - https://github.com/jsimkins2/esm_lab/blob/main/landmaskTools/ocean_mask.nc
ds = xr.open_dataset("/Users/james/Downloads/ocean_mask.nc")

ds.ocean_fraction.hvplot.quadmesh(x='lon', y='lat',project=True,geo=True, dynamic=True)

[ahuang11]

import xarray as xr #version 0.15.1
import hvplot.xarray #version 0.6.0
import holoviews as hv
import geoviews as gv #version 1.8.1
gv.extension('bokeh')

ds = xr.tutorial.open_dataset("air_temperature").isel(time=0)
da = ds["air"]

def plot(x, y):
    if x is not None:
        da_sel = da.sel(lon=x, lat=y, method="nearest")
        da.loc[dict(lat=da_sel["lat"].item(), lon=da_sel["lon"].item())] = 0
    return da.hvplot(x='lon', y='lat')

tap_stream = hv.streams.Tap()
dmap = gv.DynamicMap(plot, streams=[tap_stream])
dmap

Check out Holoviews to HoloViews — pYdeas 0.0.0 documentation for more ideas

[jsimkins2]

Thanks so much @ahuang11 ! Is there a way to get this to work with curvilinear grids? For example, when I try to use the rasm tutorial dataset, clicking removes the entire plot.

import xarray as xr #version 0.15.1
import hvplot.xarray #version 0.6.0
import holoviews as hv
import geoviews as gv #version 1.8.1
gv.extension('bokeh')

ds = xr.tutorial.open_dataset('rasm').isel(time=0)
da = ds["Tair"]

def plot(x, y):
    if x is not None:
        da_sel = da.sel(lon=x, lat=y, method="nearest")
        da.loc[dict(lat=da_sel["lat"].item(), lon=da_sel["lon"].item())] = 0
    return da.hvplot.quadmesh(x='xc', y='yc', geo=True, project=True)

tap_stream = hv.streams.Tap()
dmap = gv.DynamicMap(plot, streams=[tap_stream])
dmap

[ahuang11]

@jsimkins2

import xarray as xr #version 0.15.1
import hvplot.xarray #version 0.6.0
import holoviews as hv
import geoviews as gv #version 1.8.1
gv.extension('bokeh')

ds = xr.tutorial.open_dataset('rasm').isel(time=0)
da = ds["Tair"]

def plot(x, y):
    if x is not None:
        da_sel = da.sel(lon=x, lat=y, method="nearest")
        da.loc[dict(lat=da_sel["lat"].item(), lon=da_sel["lon"].item())] = 0
    return da.hvplot.quadmesh(x='xc', y='yc', geo=True, project=True)

tap_stream = hv.streams.Tap()
dmap = gv.DynamicMap(plot, streams=[tap_stream])
dmap

Check the browser console for errors:
“ValueError: dimensions or multi-index levels [‘lon’, ‘lat’] do not exist”

I don’t think you can use .sel with rasm since its dims are x,y not lon,lat. See how to query lat and lon in a DataArray with dims x and y? · pydata/xarray · Discussion #5009 · GitHub

[jsimkins2]

Thanks so much for your help @ahuang11 - here’s what we ended up going with for a curvilinear grid.

import xarray as xr #version 0.15.1
import hvplot.xarray #version 0.6.0
import holoviews as hv
import geoviews as gv #version 1.8.1
import cartopy.crs as ccrs
import numpy as np
gv.extension('bokeh')
clickedValues = []
# convert to numpy
# REF: https://medium.com/@petehouston/calculate-distance-of-two-locations-on-earth-using-python-1501b1944d97
def great_circle(lon1, lat1, lon2, lat2):
    #lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
    lon1 = np.radians(lon1)
    lat1 = np.radians(lat1)
    lon2 = np.radians(lon2)
    lat2 = np.radians(lat2)
    dist = 6371 * (
        np.arccos(np.sin(lat1) * np.sin(lat2) + np.cos(lat1) * np.cos(lat2) * np.cos(lon1 - lon2))
    )
    return dist

def plot(x, y):
    clickedValues.append([x,y])
    lon = x
    lat = y
    if x is not None:
        # Curvilinear
        d = great_circle(lon, lat, da['lon'], da['lat'])
        ind = np.nonzero(np.equal(d, np.amin(d)))
        yloc = int(ind[0][0])
        xloc = int(ind[1][0])
        clickedValues.append([xloc,yloc])
        da_sel = da.sel(y = yloc, x = xloc)
        clickedValues.append(dict(lat=da_sel["lat"].item(), lon=da_sel["lon"].item()))
        da.loc[yloc, xloc] = 0
    # REF: https://hvplot.holoviz.org/user_guide/Geographic_Data.html#declaring-an-output-projection
    plt = da.hvplot.quadmesh(
        'lon', 'lat', projection=ccrs.Orthographic(-90, 90),
        global_extent=True, frame_height=540, cmap='viridis',
        coastline=True
    )
    return plt

ds = xr.tutorial.open_dataset("rasm").isel(time=0)
da = ds["Tair"]
da = da.rename({
    'xc': 'lon',
    'yc': 'lat'
})
tap_stream = hv.streams.Tap()
dmap = gv.DynamicMap(plot, streams=[tap_stream])
dmap