Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multidimensional Coordinate Variables (3D lats/lons) #73

Closed
acrosby opened this issue Jun 7, 2017 · 8 comments
Closed

Multidimensional Coordinate Variables (3D lats/lons) #73

acrosby opened this issue Jun 7, 2017 · 8 comments

Comments

@acrosby
Copy link

acrosby commented Jun 7, 2017

We have been working on CF-compliant netcdf formats to provide flexible multi resolution (but gridded) input to unstructured hydrodynamic and wave models. xarray appears to support dealing with these files/groups relatively painlessly, and it would be great if we could get support in geoviews which would allow us to demonstrate the capabilities of the format nicely.

Here is an example of a group within a test file:

<xarray.Dataset>
Dimensions:  (time: 168, xi: 25, yi: 25)
Coordinates:
  * time     (time) datetime64[ns] 2017-05-31T12:00:00 2017-05-31T13:00:00 ...
    lon      (time, yi, xi) float64 -180.0 -165.0 -150.0 -135.0 -120.0 ...
    lat      (time, yi, xi) float64 -70.0 -70.0 -70.0 -70.0 -70.0 -70.0 ...
Dimensions without coordinates: xi, yi
Data variables:
    U10      (time, yi, xi) float64 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 ...
    V10      (time, yi, xi) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
    PSFC     (time, yi, xi) float64 1.015e+03 1.015e+03 1.015e+03 1.015e+03 ...
Attributes:
    rank: 1

This is the test I was attempting:

xr = gv.Dataset(xarray.open_dataset("./test_file.nc", group="region"), 
                kdims=["lon", "lat", "time"], vdims=["PSFC", "U10", "V10"], 
                crs=ccrs.PlateCarree())
xr.to(gv.Image, ['lon', 'lat'], "PSFC", dynamic=True) * gf.coastline

Here is the full traceback (although I understand from Gitter that there are probably bigger issues at play to support a structure like this):

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/home/alex/.local/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
    305                 pass
    306             else:
--> 307                 return printer(obj)
    308             # Finally look for special method names
    309             method = get_real_method(obj, self.print_method)

/usr/local/lib/python2.7/dist-packages/holoviews/ipython/display_hooks.pyc in pprint_display(obj)
    255     if not ip.display_formatter.formatters['text/plain'].pprint:
    256         return None
--> 257     return display(obj, raw=True)
    258 
    259 

/usr/local/lib/python2.7/dist-packages/holoviews/ipython/display_hooks.pyc in display(obj, raw, **kwargs)
    241     elif isinstance(obj, (HoloMap, DynamicMap)):
    242         with option_state(obj):
--> 243             html = map_display(obj)
    244     else:
    245         return repr(obj) if raw else IPython.display.display(obj, **kwargs)

/usr/local/lib/python2.7/dist-packages/holoviews/ipython/display_hooks.pyc in wrapped(element)
    127         try:
    128             html = fn(element,
--> 129                       max_frames=OutputMagic.options['max_frames'])
    130 
    131             # Only want to add to the archive for one display hook...

/usr/local/lib/python2.7/dist-packages/holoviews/ipython/display_hooks.pyc in map_display(vmap, max_frames)
    196         return None
    197 
--> 198     return render(vmap)
    199 
    200 

/usr/local/lib/python2.7/dist-packages/holoviews/ipython/display_hooks.pyc in render(obj, **kwargs)
     57     if renderer.fig == 'pdf':
     58         renderer = renderer.instance(fig='png')
---> 59     return renderer.html(obj, **kwargs)
     60 
     61 

/usr/local/lib/python2.7/dist-packages/holoviews/plotting/renderer.pyc in html(self, obj, fmt, css, comm, **kwargs)
    253         code to initialize a Comm, if the plot supplies one.
    254         """
--> 255         plot, fmt =  self._validate(obj, fmt)
    256         figdata, _ = self(plot, fmt, **kwargs)
    257         if css is None: css = self.css

/usr/local/lib/python2.7/dist-packages/holoviews/plotting/renderer.pyc in _validate(self, obj, fmt)
    189         if isinstance(obj, tuple(self.widgets.values())):
    190             return obj, 'html'
--> 191         plot = self.get_plot(obj, renderer=self)
    192 
    193         fig_formats = self.mode_formats['fig'][self.mode]

/usr/local/lib/python2.7/dist-packages/holoviews/plotting/renderer.pyc in get_plot(self_or_cls, obj, renderer)
    164         """
    165         # Initialize DynamicMaps with first data item
--> 166         initialize_dynamic(obj)
    167 
    168         if not isinstance(obj, Plot) and not displayable(obj):

/usr/local/lib/python2.7/dist-packages/holoviews/plotting/util.pyc in initialize_dynamic(obj)
    173             continue
    174         if not len(dmap):
--> 175             dmap[dmap._initial_key()]
    176 
    177 

/usr/local/lib/python2.7/dist-packages/holoviews/core/spaces.pyc in __getitem__(self, key)
    942         # Not a cross product and nothing cached so compute element.
    943         if cache is not None: return cache
--> 944         val = self._execute_callback(*tuple_key)
    945         if data_slice:
    946             val = self._dataslice(val, data_slice)

/usr/local/lib/python2.7/dist-packages/holoviews/core/spaces.pyc in _execute_callback(self, *args)
    791 
    792         with dynamicmap_memoization(self.callback, self.streams):
--> 793             retval = self.callback(*args, **kwargs)
    794         return self._style(retval)
    795 

/usr/local/lib/python2.7/dist-packages/holoviews/core/spaces.pyc in __call__(self, *args, **kwargs)
    519 
    520         try:
--> 521             ret = self.callable(*args, **kwargs)
    522         except:
    523             posstr = ', '.join(['%r' % el for el in inargs]) if inargs else ''

/usr/local/lib/python2.7/dist-packages/holoviews/core/spaces.pyc in dynamic_mul(*args, **kwargs)
    207             if isinstance(self, DynamicMap):
    208                 def dynamic_mul(*args, **kwargs):
--> 209                     element = self[args]
    210                     return element * other
    211                 callback = Callable(dynamic_mul, inputs=[self, other])

/usr/local/lib/python2.7/dist-packages/holoviews/core/spaces.pyc in __getitem__(self, key)
    942         # Not a cross product and nothing cached so compute element.
    943         if cache is not None: return cache
--> 944         val = self._execute_callback(*tuple_key)
    945         if data_slice:
    946             val = self._dataslice(val, data_slice)

/usr/local/lib/python2.7/dist-packages/holoviews/core/spaces.pyc in _execute_callback(self, *args)
    791 
    792         with dynamicmap_memoization(self.callback, self.streams):
--> 793             retval = self.callback(*args, **kwargs)
    794         return self._style(retval)
    795 

/usr/local/lib/python2.7/dist-packages/holoviews/core/spaces.pyc in __call__(self, *args, **kwargs)
    519 
    520         try:
--> 521             ret = self.callable(*args, **kwargs)
    522         except:
    523             posstr = ', '.join(['%r' % el for el in inargs]) if inargs else ''

/usr/local/lib/python2.7/dist-packages/holoviews/core/data/__init__.pyc in load_subset(*args)
    543                 if drop_dim and self.interface.gridded:
    544                     data = data.columns()
--> 545                 return group_type(data, **group_kwargs)
    546             dynamic_dims = [d(values=list(self.interface.values(self, d.name, False)))
    547                             for d in dimensions]

/usr/local/lib/python2.7/dist-packages/geoviews/element/geo.pyc in __init__(self, data, **kwargs)
     81         elif crs:
     82             kwargs['crs'] = crs
---> 83         super(_Element, self).__init__(data, **kwargs)
     84 
     85 

/usr/local/lib/python2.7/dist-packages/holoviews/element/raster.pyc in __init__(self, data, bounds, extents, xdensity, ydensity, **params)
    242         if bounds is None:
    243             xvals = self.dimension_values(0, False)
--> 244             l, r, xdensity, _ = util.bound_range(xvals, xdensity)
    245             yvals = self.dimension_values(1, False)
    246             b, t, ydensity, _ = util.bound_range(yvals, ydensity)

/usr/local/lib/python2.7/dist-packages/holoviews/core/util.pyc in bound_range(vals, density)
   1375     low, high = vals.min(), vals.max()
   1376     invert = False
-> 1377     if vals[0] > vals[1]:
   1378         invert = True
   1379     if not density:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

:DynamicMap   [time]

Versions I am currently using if it's relevant:

geoviews==1.2.0
holoviews==1.7.0
xarray==0.9.1
@philippjfr
Copy link
Member

Should have realized this earlier but this is related to:

#57

It's possible that we can get the interfaces to interpret 2d coordinate arrays correctly. The main issue is in plotting this data in bokeh since the data will have to be regridded before being rendered with bokeh.

@philippjfr
Copy link
Member

I've put together a notebook outlining how to work with datasets like this in GeoViews:

https://anaconda.org/philippjfr/working_with_multi-dimensional_coordinates/notebook

@snowman2
Copy link

snowman2 commented Oct 5, 2017

@philippjfr, the example you gave already starts with 1D x,y coordinates. Do you have an example for a dataset with only 2D coordinates? It would be nice to see a HoloViews version of this example: http://xarray.pydata.org/en/stable/examples/multidimensional-coords.html#plotting

@philippjfr
Copy link
Member

philippjfr commented Oct 5, 2017

I'll try to put together an example. That said in the long term we should implement an actual interface for 2D coordinates, so you could easily do something like this:

ds = hv.Dataset(xr.tutorial.load_dataset('rasm'), ['x', 'y', 'time'], 'Tair')
datashade(ds.groupby('time'))

@snowman2
Copy link

snowman2 commented Oct 5, 2017

<xarray.Dataset>
Dimensions:                                         (time: 1, x: 3856, y: 1624)
Coordinates:
  * time                                            (time) datetime64[ns] 2017-10-02T23:58:55.816000 ...
    cell_lat                                        (y, x) float64 84.66 ...
    cell_lon                                        (y, x) float64 -180.0 ...
Dimensions without coordinates: x, y
Data variables:
    Forecast_Data/tb_h_forecast_ensstd              (y, x) float64 nan nan ...
    Analysis_Data/surface_temp_analysis             (y, x) float64 nan nan ...
    Observations_Data/tb_h_obs_errstd               (y, x) float64 nan nan ...
    Observations_Data/tb_v_obs_errstd               (y, x) float64 nan nan ...
    Forecast_Data/surface_temp_forecast             (y, x) float64 nan nan ...
    Analysis_Data/sm_surface_analysis               (y, x) float64 nan nan ...
    Forecast_Data/sm_surface_forecast              (y, x) float64 nan nan ...
%%opts Overlay [width=600 height=400] Image (cmap='viridis') Feature (line_color='black')
sm_surf = gv.Dataset(sa_xds, 
                     kdims=['cell_lon', 'cell_lat'],
                     group='Forecast_Data/sm_surface_forecast', 
                     vdims=['Forecast_Data/sm_surface_forecast'],
                     crs=ccrs.PlateCarree())
sm_img = sm_surf.to(gv.Image, ['cell_lon', 'cell_lat'], 'Forecast_Data/sm_surface_forecast', dynamic=True)
WARNING:root:Image07259: Setting non-parameter attribute dynamic=True using a mechanism intended only for parameters

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-22-7b037eec2c19> in <module>()
      4                      vdims=['Forecast_Data/sm_surface_forecast'],
      5                      crs=ccrs.PlateCarree())
----> 6 sm_img = sm_surf.to(gv.Image, ['cell_lon', 'cell_lat'], 'Forecast_Data/sm_surface_forecast', dynamic=True)

/miniconda3/envs/dp2/lib/python2.7/site-packages/geoviews/element/__init__.pyc in __call__(self, *args, **kwargs)
     24         if 'crs' not in kwargs and issubclass(group_type, _Element):
     25             kwargs['crs'] = self._element.crs
---> 26         return super(GeoConversion, self).__call__(*args, **kwargs)
     27 
     28     def linecontours(self, kdims=None, vdims=None, mdims=None, **kwargs):

/miniconda3/envs/dp2/lib/python2.7/site-packages/holoviews/core/data/__init__.pyc in __call__(self, new_type, kdims, vdims, groupby, sort, **kwargs)
    129         params.update(kwargs)
    130         if len(kdims) == selected.ndims or not groupby:
--> 131             element = new_type(selected, **params)
    132             return element.sort() if sort else element
    133         group = selected.groupby(groupby, container_type=HoloMap,

/miniconda3/envs/dp2/lib/python2.7/site-packages/geoviews/element/geo.pyc in __init__(self, data, **kwargs)
     87         elif crs:
     88             kwargs['crs'] = crs
---> 89         super(_Element, self).__init__(data, **kwargs)
     90 
     91 

/miniconda3/envs/dp2/lib/python2.7/site-packages/holoviews/element/raster.pyc in __init__(self, data, bounds, extents, xdensity, ydensity, **params)
    242         if bounds is None:
    243             xvals = self.dimension_values(0, False)
--> 244             l, r, xdensity, _ = util.bound_range(xvals, xdensity)
    245             yvals = self.dimension_values(1, False)
    246             b, t, ydensity, _ = util.bound_range(yvals, ydensity)

/miniconda3/envs/dp2/lib/python2.7/site-packages/holoviews/core/util.pyc in bound_range(vals, density)
   1480     low, high = vals.min(), vals.max()
   1481     invert = False
-> 1482     if len(vals) > 1 and vals[0] > vals[1]:
   1483         invert = True
   1484     if not density:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
Name: geoviews
Version: 1.3.2
Name: holoviews
Version: 1.8.4

@snowman2
Copy link

snowman2 commented Oct 5, 2017

This is an example of a solution that may help with the error:
https://github.com/snowman2/pangaea/blob/master/pangaea/xlsm.py#L91-L102

@philippjfr
Copy link
Member

philippjfr commented Jan 6, 2018

Getting somewhere here, once #116 is merged you'll be able to do something like this:

%%opts Image [width=700 height=350 clipping_colors={'NaN': (0,0,0,0)} colorbar=True]
tiles = gv.WMTS('https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{Z}/{Y}/{X}.jpg')

xrds = xr.tutorial.load_dataset('RASM_example_data')
qmesh = gv.Dataset(xrds, ['time', 'xc', 'yc']).to(gv.QuadMesh, ['xc', 'yc'])

tiles * rasterize(gv.operation.project(qmesh), aggregator=ds.mean('Tair'),
                  precompute=True)

qmesh_regrid

@philippjfr
Copy link
Member

This is now merged, closing.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants