Resampling Grids

import numpy as np
import xarray as xr
import geoviews as gv
import datashader as dsh
from geoviews import opts

gv.extension('bokeh', 'matplotlib')

opts.defaults(
    opts.Image(width=600, height=400, colorbar=True),
    opts.Feature(apply_ranges=False),
    opts.QuadMesh(width=600, height=400, colorbar=True))

In geographical applications grids and meshes of different kinds are very common and for visualization and analysis it is often very important to be able to resample them in different ways. Regridding can refer both to upsampling and downsampling a grid or mesh, which is achieved through interpolation and aggregation.

Naive approaches to regridding treat the space as flat, which is often simpler but can also give less accurate results when working with a spherical space such as the earth. In this user guide we will summarize how to work with different grid types including rectilinear, curvilinear grids and trimeshes. Additionally we will discuss different approaches to regridding working based on the assumption of a flat earth (using datashader) and a spherical earth (xESMF).

Rectilinear grids

Rectilinear grids are one of the most standard formats and are defined by regularly sampled coordinates along the two axes. The air_temperature dataset provided by xarray and used throughout the GeoViews documentation provides a good example.

ds = xr.tutorial.open_dataset('air_temperature').load().isel(time=slice(0, 100))
ds

<xarray.Dataset> Size: 1MB
Dimensions:  (lat: 25, time: 100, lon: 53)
Coordinates:
  * lat      (lat) float32 100B 75.0 72.5 70.0 67.5 65.0 ... 22.5 20.0 17.5 15.0
  * lon      (lon) float32 212B 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0
  * time     (time) datetime64[ns] 800B 2013-01-01 ... 2013-01-25T18:00:00
Data variables:
    air      (time, lat, lon) float64 1MB 241.2 242.5 243.5 ... 295.4 294.9
Attributes:
    Conventions:  COARDS
    title:        4x daily NMC reanalysis (1948)
    description:  Data is from NMC initialized reanalysis\n(4x/day).  These a...
    platform:     Model
    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...

xarray.Dataset

• Dimensions:
  • lat: 25
  • time: 100
  • lon: 53
• Coordinates: (3)
  • lat
    (lat)
    float32
    75.0 72.5 70.0 ... 20.0 17.5 15.0

    standard_name :

    latitude

    long_name :

    Latitude

    units :

    degrees_north

    axis :

    Y

    array([75. , 72.5, 70. , 67.5, 65. , 62.5, 60. , 57.5, 55. , 52.5, 50. , 47.5,
           45. , 42.5, 40. , 37.5, 35. , 32.5, 30. , 27.5, 25. , 22.5, 20. , 17.5,
           15. ], dtype=float32)

  • lon
    (lon)
    float32
    200.0 202.5 205.0 ... 327.5 330.0

    standard_name :

    longitude

    long_name :

    Longitude

    units :

    degrees_east

    axis :

    X

    array([200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
           225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
           250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
           275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
           300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
           325. , 327.5, 330. ], dtype=float32)

  • time
    (time)
    datetime64[ns]
    2013-01-01 ... 2013-01-25T18:00:00

    standard_name :

    time

    long_name :

    Time

    array(['2013-01-01T00:00:00.000000000', '2013-01-01T06:00:00.000000000',
           '2013-01-01T12:00:00.000000000', '2013-01-01T18:00:00.000000000',
           '2013-01-02T00:00:00.000000000', '2013-01-02T06:00:00.000000000',
           '2013-01-02T12:00:00.000000000', '2013-01-02T18:00:00.000000000',
           '2013-01-03T00:00:00.000000000', '2013-01-03T06:00:00.000000000',
           '2013-01-03T12:00:00.000000000', '2013-01-03T18:00:00.000000000',
           '2013-01-04T00:00:00.000000000', '2013-01-04T06:00:00.000000000',
           '2013-01-04T12:00:00.000000000', '2013-01-04T18:00:00.000000000',
           '2013-01-05T00:00:00.000000000', '2013-01-05T06:00:00.000000000',
           '2013-01-05T12:00:00.000000000', '2013-01-05T18:00:00.000000000',
           '2013-01-06T00:00:00.000000000', '2013-01-06T06:00:00.000000000',
           '2013-01-06T12:00:00.000000000', '2013-01-06T18:00:00.000000000',
           '2013-01-07T00:00:00.000000000', '2013-01-07T06:00:00.000000000',
           '2013-01-07T12:00:00.000000000', '2013-01-07T18:00:00.000000000',
           '2013-01-08T00:00:00.000000000', '2013-01-08T06:00:00.000000000',
           '2013-01-08T12:00:00.000000000', '2013-01-08T18:00:00.000000000',
           '2013-01-09T00:00:00.000000000', '2013-01-09T06:00:00.000000000',
           '2013-01-09T12:00:00.000000000', '2013-01-09T18:00:00.000000000',
           '2013-01-10T00:00:00.000000000', '2013-01-10T06:00:00.000000000',
           '2013-01-10T12:00:00.000000000', '2013-01-10T18:00:00.000000000',
           '2013-01-11T00:00:00.000000000', '2013-01-11T06:00:00.000000000',
           '2013-01-11T12:00:00.000000000', '2013-01-11T18:00:00.000000000',
           '2013-01-12T00:00:00.000000000', '2013-01-12T06:00:00.000000000',
           '2013-01-12T12:00:00.000000000', '2013-01-12T18:00:00.000000000',
           '2013-01-13T00:00:00.000000000', '2013-01-13T06:00:00.000000000',
           '2013-01-13T12:00:00.000000000', '2013-01-13T18:00:00.000000000',
           '2013-01-14T00:00:00.000000000', '2013-01-14T06:00:00.000000000',
           '2013-01-14T12:00:00.000000000', '2013-01-14T18:00:00.000000000',
           '2013-01-15T00:00:00.000000000', '2013-01-15T06:00:00.000000000',
           '2013-01-15T12:00:00.000000000', '2013-01-15T18:00:00.000000000',
           '2013-01-16T00:00:00.000000000', '2013-01-16T06:00:00.000000000',
           '2013-01-16T12:00:00.000000000', '2013-01-16T18:00:00.000000000',
           '2013-01-17T00:00:00.000000000', '2013-01-17T06:00:00.000000000',
           '2013-01-17T12:00:00.000000000', '2013-01-17T18:00:00.000000000',
           '2013-01-18T00:00:00.000000000', '2013-01-18T06:00:00.000000000',
           '2013-01-18T12:00:00.000000000', '2013-01-18T18:00:00.000000000',
           '2013-01-19T00:00:00.000000000', '2013-01-19T06:00:00.000000000',
           '2013-01-19T12:00:00.000000000', '2013-01-19T18:00:00.000000000',
           '2013-01-20T00:00:00.000000000', '2013-01-20T06:00:00.000000000',
           '2013-01-20T12:00:00.000000000', '2013-01-20T18:00:00.000000000',
           '2013-01-21T00:00:00.000000000', '2013-01-21T06:00:00.000000000',
           '2013-01-21T12:00:00.000000000', '2013-01-21T18:00:00.000000000',
           '2013-01-22T00:00:00.000000000', '2013-01-22T06:00:00.000000000',
           '2013-01-22T12:00:00.000000000', '2013-01-22T18:00:00.000000000',
           '2013-01-23T00:00:00.000000000', '2013-01-23T06:00:00.000000000',
           '2013-01-23T12:00:00.000000000', '2013-01-23T18:00:00.000000000',
           '2013-01-24T00:00:00.000000000', '2013-01-24T06:00:00.000000000',
           '2013-01-24T12:00:00.000000000', '2013-01-24T18:00:00.000000000',
           '2013-01-25T00:00:00.000000000', '2013-01-25T06:00:00.000000000',
           '2013-01-25T12:00:00.000000000', '2013-01-25T18:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (1)
  • air
    (time, lat, lon)
    float64
    241.2 242.5 243.5 ... 295.4 294.9

    long_name :

    4xDaily Air temperature at sigma level 995

    units :

    degK

    precision :

    2

    GRIB_id :

    11

    GRIB_name :

    TMP

    var_desc :

    Air temperature

    dataset :

    NMC Reanalysis

    level_desc :

    Surface

    statistic :

    Individual Obs

    parent_stat :

    Other

    actual_range :

    [185.16 322.1 ]

    array([[[241.2 , 242.5 , 243.5 , ..., 232.8 , 235.5 , 238.6 ],
            [243.8 , 244.5 , 244.7 , ..., 232.8 , 235.3 , 239.3 ],
            [250.  , 249.8 , 248.89, ..., 233.2 , 236.39, 241.7 ],
            ...,
            [296.6 , 296.2 , 296.4 , ..., 295.4 , 295.1 , 294.7 ],
            [295.9 , 296.2 , 296.79, ..., 295.9 , 295.9 , 295.2 ],
            [296.29, 296.79, 297.1 , ..., 296.9 , 296.79, 296.6 ]],
    
           [[242.1 , 242.7 , 243.1 , ..., 232.  , 233.6 , 235.8 ],
            [243.6 , 244.1 , 244.2 , ..., 231.  , 232.5 , 235.7 ],
            [253.2 , 252.89, 252.1 , ..., 230.8 , 233.39, 238.5 ],
            ...,
            [296.4 , 295.9 , 296.2 , ..., 295.4 , 295.1 , 294.79],
            [296.2 , 296.7 , 296.79, ..., 295.6 , 295.5 , 295.1 ],
            [296.29, 297.2 , 297.4 , ..., 296.4 , 296.4 , 296.6 ]],
    
           [[242.3 , 242.2 , 242.3 , ..., 234.3 , 236.1 , 238.7 ],
            [244.6 , 244.39, 244.  , ..., 230.3 , 232.  , 235.7 ],
            [256.2 , 255.5 , 254.2 , ..., 231.2 , 233.2 , 238.2 ],
            ...,
    ...
            ...,
            [296.  , 295.4 , 295.2 , ..., 294.  , 293.4 , 292.4 ],
            [297.  , 296.7 , 296.1 , ..., 293.9 , 293.6 , 292.9 ],
            [297.29, 297.29, 297.  , ..., 294.7 , 294.6 , 294.5 ]],
    
           [[251.8 , 252.89, 253.5 , ..., 242.  , 242.  , 242.5 ],
            [254.2 , 254.3 , 253.7 , ..., 241.2 , 241.  , 241.89],
            [248.8 , 248.  , 246.7 , ..., 239.7 , 240.89, 243.5 ],
            ...,
            [295.  , 294.  , 293.7 , ..., 294.2 , 293.5 , 292.7 ],
            [296.2 , 295.7 , 295.1 , ..., 294.2 , 294.2 , 293.29],
            [296.79, 296.79, 296.5 , ..., 295.1 , 295.  , 294.7 ]],
    
           [[249.3 , 250.5 , 251.2 , ..., 240.6 , 240.7 , 241.3 ],
            [252.89, 254.1 , 254.5 , ..., 240.39, 240.3 , 241.3 ],
            [250.7 , 250.6 , 250.  , ..., 238.7 , 240.  , 242.6 ],
            ...,
            [295.9 , 294.5 , 294.5 , ..., 294.2 , 293.1 , 292.7 ],
            [296.1 , 295.79, 295.6 , ..., 294.6 , 294.29, 293.9 ],
            [296.79, 296.79, 296.6 , ..., 295.79, 295.4 , 294.9 ]]])

• Indexes: (3)
  • lat
    PandasIndex

    PandasIndex(Index([75.0, 72.5, 70.0, 67.5, 65.0, 62.5, 60.0, 57.5, 55.0, 52.5, 50.0, 47.5,
           45.0, 42.5, 40.0, 37.5, 35.0, 32.5, 30.0, 27.5, 25.0, 22.5, 20.0, 17.5,
           15.0],
          dtype='float32', name='lat'))

  • lon
    PandasIndex

    PandasIndex(Index([200.0, 202.5, 205.0, 207.5, 210.0, 212.5, 215.0, 217.5, 220.0, 222.5,
           225.0, 227.5, 230.0, 232.5, 235.0, 237.5, 240.0, 242.5, 245.0, 247.5,
           250.0, 252.5, 255.0, 257.5, 260.0, 262.5, 265.0, 267.5, 270.0, 272.5,
           275.0, 277.5, 280.0, 282.5, 285.0, 287.5, 290.0, 292.5, 295.0, 297.5,
           300.0, 302.5, 305.0, 307.5, 310.0, 312.5, 315.0, 317.5, 320.0, 322.5,
           325.0, 327.5, 330.0],
          dtype='float32', name='lon'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2013-01-01 00:00:00', '2013-01-01 06:00:00',
                   '2013-01-01 12:00:00', '2013-01-01 18:00:00',
                   '2013-01-02 00:00:00', '2013-01-02 06:00:00',
                   '2013-01-02 12:00:00', '2013-01-02 18:00:00',
                   '2013-01-03 00:00:00', '2013-01-03 06:00:00',
                   '2013-01-03 12:00:00', '2013-01-03 18:00:00',
                   '2013-01-04 00:00:00', '2013-01-04 06:00:00',
                   '2013-01-04 12:00:00', '2013-01-04 18:00:00',
                   '2013-01-05 00:00:00', '2013-01-05 06:00:00',
                   '2013-01-05 12:00:00', '2013-01-05 18:00:00',
                   '2013-01-06 00:00:00', '2013-01-06 06:00:00',
                   '2013-01-06 12:00:00', '2013-01-06 18:00:00',
                   '2013-01-07 00:00:00', '2013-01-07 06:00:00',
                   '2013-01-07 12:00:00', '2013-01-07 18:00:00',
                   '2013-01-08 00:00:00', '2013-01-08 06:00:00',
                   '2013-01-08 12:00:00', '2013-01-08 18:00:00',
                   '2013-01-09 00:00:00', '2013-01-09 06:00:00',
                   '2013-01-09 12:00:00', '2013-01-09 18:00:00',
                   '2013-01-10 00:00:00', '2013-01-10 06:00:00',
                   '2013-01-10 12:00:00', '2013-01-10 18:00:00',
                   '2013-01-11 00:00:00', '2013-01-11 06:00:00',
                   '2013-01-11 12:00:00', '2013-01-11 18:00:00',
                   '2013-01-12 00:00:00', '2013-01-12 06:00:00',
                   '2013-01-12 12:00:00', '2013-01-12 18:00:00',
                   '2013-01-13 00:00:00', '2013-01-13 06:00:00',
                   '2013-01-13 12:00:00', '2013-01-13 18:00:00',
                   '2013-01-14 00:00:00', '2013-01-14 06:00:00',
                   '2013-01-14 12:00:00', '2013-01-14 18:00:00',
                   '2013-01-15 00:00:00', '2013-01-15 06:00:00',
                   '2013-01-15 12:00:00', '2013-01-15 18:00:00',
                   '2013-01-16 00:00:00', '2013-01-16 06:00:00',
                   '2013-01-16 12:00:00', '2013-01-16 18:00:00',
                   '2013-01-17 00:00:00', '2013-01-17 06:00:00',
                   '2013-01-17 12:00:00', '2013-01-17 18:00:00',
                   '2013-01-18 00:00:00', '2013-01-18 06:00:00',
                   '2013-01-18 12:00:00', '2013-01-18 18:00:00',
                   '2013-01-19 00:00:00', '2013-01-19 06:00:00',
                   '2013-01-19 12:00:00', '2013-01-19 18:00:00',
                   '2013-01-20 00:00:00', '2013-01-20 06:00:00',
                   '2013-01-20 12:00:00', '2013-01-20 18:00:00',
                   '2013-01-21 00:00:00', '2013-01-21 06:00:00',
                   '2013-01-21 12:00:00', '2013-01-21 18:00:00',
                   '2013-01-22 00:00:00', '2013-01-22 06:00:00',
                   '2013-01-22 12:00:00', '2013-01-22 18:00:00',
                   '2013-01-23 00:00:00', '2013-01-23 06:00:00',
                   '2013-01-23 12:00:00', '2013-01-23 18:00:00',
                   '2013-01-24 00:00:00', '2013-01-24 06:00:00',
                   '2013-01-24 12:00:00', '2013-01-24 18:00:00',
                   '2013-01-25 00:00:00', '2013-01-25 06:00:00',
                   '2013-01-25 12:00:00', '2013-01-25 18:00:00'],
                  dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (5)

  Conventions :

  COARDS

  title :

  4x daily NMC reanalysis (1948)

  description :

  Data is from NMC initialized reanalysis (4x/day). These are the 0.9950 sigma level values.

  platform :

  Model

  references :

  http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html

As we already know from the Gridded Dataset sections, an xarray of this kind can easily be wrapped in a GeoViews dataset:

gvds = gv.Dataset(ds).redim.range(air=(230, 300))
gvds

:Dataset   [lat,lon,time]   (4xDaily Air temperature at sigma level 995)

We can also easily plot this data by using the .to method to group the data into a set of Image elements indexed by ‘time’:

images = gvds.to(gv.Image, ['lon', 'lat'], dynamic=True)

It is important to note that if we look at the longitude coordinates above we can see that they are defined in the range (0, 360), while GeoViews generally expects it to be in the range (-180, 180). To correct this we can apply a simple fix:

ds['lon'] = np.where(ds.lon>180, ds.lon-360, ds.lon)

Now we are ready to display the data. Note that throughout this user guide we will be using Bokeh but we could easily switch to matplotlib if needed.

opts.defaults(opts.Image(cmap='viridis'))

images * gv.feature.coastline

Datashader

HoloViews provides high-level wrappers around the datashader library, which make it possible to quickly resample or aggregate datasets of different kinds. Datashader knows nothing about non-flat coordinate systems, but provides a very fast, parallelized regridding operation for rectilinear grids. Here we will import the regrid operation and pass it our stack of images from above. While this dataset is fairly small and regridding will actually upsample the image to match the dimensions of the plot, regrid can very quickly downsample very large datasets.

One important thing to note about the resampling operations we will be working with in this user guide is that they are dynamic and linked to the plot dimensions and axis ranges. This means that whenever we zoom or pan the data will be resampled. If we want to disable this linked behavior and supply an explicit width and height we can disable the streams by passing streams=[] as a keyword argument.

from holoviews.operation.datashader import regrid
regrid(images) * gv.feature.coastline

xESMF

The xESMF library is specifically designed to provide an easy way to accurately resample grids defined in geographic coordinate systems and differs significantly from the simpler approach used by datashader, which applies simple upsampling and downsampling. xESMF is a wrapper around the ESMF regridding algorithms, which compute an interpolation weight matrix which is applied to remap the values of the source grid onto the destination grid.

In GeoViews these algorithms are made available via the weighted_regrid operation, which supports the different interpolation modes including: ‘bilinear’, ‘nearest_s2d’, ‘nearest_d2s’ and ‘conservative’. Since generating the sparse weight matrix takes much longer than applying it the operation will cache the weight matrix on disk for later use; this optimization can be disabled via the reuse_weights parameter or customized by defining a custom file_pattern.

from geoviews.operation.regrid import weighted_regrid
weighted_regrid(images) * gv.feature.coastline

Since this operation creates local weight files we will want to clean up after ourselves once we are done, to do so we can call the weighted_regrid.clean_weight_files method.

weighted_regrid.clean_weight_files()

Deleted 0 weight files

Onto an existing grid

The operation also allows us to define a target grid, which we can either define manually or by using a utility provided by the xESMF library. Here we will define a $2^{\circ}\times2^{\circ}$ grid.

import xesmf as xe
grid = xe.util.grid_2d(-160, -35, 2, 15, 70, 2)
grid

<xarray.Dataset> Size: 55kB
Dimensions:  (y: 27, x: 62, y_b: 28, x_b: 63)
Coordinates:
    lon      (y, x) float64 13kB -159.0 -157.0 -155.0 ... -41.0 -39.0 -37.0
    lat      (y, x) float64 13kB 16.0 16.0 16.0 16.0 ... 68.0 68.0 68.0 68.0
    lon_b    (y_b, x_b) float64 14kB -160.0 -158.0 -156.0 ... -40.0 -38.0 -36.0
    lat_b    (y_b, x_b) float64 14kB 15.0 15.0 15.0 15.0 ... 69.0 69.0 69.0 69.0
Dimensions without coordinates: y, x, y_b, x_b
Data variables:
    *empty*

xarray.Dataset

• Dimensions:
  • y: 27
  • x: 62
  • y_b: 28
  • x_b: 63
• Coordinates: (4)
  • lon
    (y, x)
    float64
    -159.0 -157.0 ... -39.0 -37.0

    standard_name :

    longitude

    array([[-159., -157., -155., ...,  -41.,  -39.,  -37.],
           [-159., -157., -155., ...,  -41.,  -39.,  -37.],
           [-159., -157., -155., ...,  -41.,  -39.,  -37.],
           ...,
           [-159., -157., -155., ...,  -41.,  -39.,  -37.],
           [-159., -157., -155., ...,  -41.,  -39.,  -37.],
           [-159., -157., -155., ...,  -41.,  -39.,  -37.]])

  • lat
    (y, x)
    float64
    16.0 16.0 16.0 ... 68.0 68.0 68.0

    standard_name :

    latitude

    array([[16., 16., 16., ..., 16., 16., 16.],
           [18., 18., 18., ..., 18., 18., 18.],
           [20., 20., 20., ..., 20., 20., 20.],
           ...,
           [64., 64., 64., ..., 64., 64., 64.],
           [66., 66., 66., ..., 66., 66., 66.],
           [68., 68., 68., ..., 68., 68., 68.]])

  • lon_b
    (y_b, x_b)
    float64
    -160.0 -158.0 ... -38.0 -36.0

    array([[-160., -158., -156., ...,  -40.,  -38.,  -36.],
           [-160., -158., -156., ...,  -40.,  -38.,  -36.],
           [-160., -158., -156., ...,  -40.,  -38.,  -36.],
           ...,
           [-160., -158., -156., ...,  -40.,  -38.,  -36.],
           [-160., -158., -156., ...,  -40.,  -38.,  -36.],
           [-160., -158., -156., ...,  -40.,  -38.,  -36.]])

  • lat_b
    (y_b, x_b)
    float64
    15.0 15.0 15.0 ... 69.0 69.0 69.0

    array([[15., 15., 15., ..., 15., 15., 15.],
           [17., 17., 17., ..., 17., 17., 17.],
           [19., 19., 19., ..., 19., 19., 19.],
           ...,
           [65., 65., 65., ..., 65., 65., 65.],
           [67., 67., 67., ..., 67., 67., 67.],
           [69., 69., 69., ..., 69., 69., 69.]])

• Data variables: (0)
• Indexes: (0)
• Attributes: (0)

Since the grid has 2D coordinate arrays the regridded data will be wrapped in and displayed as a QuadMesh:

target = gv.Dataset(grid, kdims=['lon', 'lat'])
weighted_regrid(images, target=target, streams=[]) * gv.feature.coastline

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

Curvilinear Grids

Curvilinear grids are another very common mesh type, which are usually defined by multi-dimensional coordinate arrays:

ds = xr.tutorial.open_dataset('rasm').load()
ds

<xarray.Dataset> Size: 17MB
Dimensions:  (time: 36, y: 205, x: 275)
Coordinates:
  * time     (time) object 288B 1980-09-16 12:00:00 ... 1983-08-17 00:00:00
    xc       (y, x) float64 451kB 189.2 189.4 189.6 189.7 ... 17.4 17.15 16.91
    yc       (y, x) float64 451kB 16.53 16.78 17.02 17.27 ... 28.01 27.76 27.51
Dimensions without coordinates: y, x
Data variables:
    Tair     (time, y, x) float64 16MB nan nan nan nan ... 28.66 28.19 28.21
Attributes:
    title:                     /workspace/jhamman/processed/R1002RBRxaaa01a/l...
    institution:               U.W.
    source:                    RACM R1002RBRxaaa01a
    output_frequency:          daily
    output_mode:               averaged
    convention:                CF-1.4
    references:                Based on the initial model of Liang et al., 19...
    comment:                   Output from the Variable Infiltration Capacity...
    nco_openmp_thread_number:  1
    NCO:                       netCDF Operators version 4.7.9 (Homepage = htt...
    history:                   Fri Aug  7 17:57:38 2020: ncatted -a bounds,,d...

xarray.Dataset

• Dimensions:
  • time: 36
  • y: 205
  • x: 275
• Coordinates: (3)
  • time
    (time)
    object
    1980-09-16 12:00:00 ... 1983-08-...

    long_name :

    time

    type_preferred :

    int

    array([cftime.DatetimeNoLeap(1980, 9, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1980, 10, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1980, 11, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1980, 12, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 1, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 2, 15, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 3, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 4, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 5, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 6, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 7, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 8, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 9, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 10, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 11, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1981, 12, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 1, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 2, 15, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 3, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 4, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 5, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 6, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 7, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 8, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 9, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 10, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 11, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1982, 12, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1983, 1, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1983, 2, 15, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1983, 3, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1983, 4, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1983, 5, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1983, 6, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1983, 7, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1983, 8, 17, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • xc
    (y, x)
    float64
    189.2 189.4 189.6 ... 17.15 16.91

    long_name :

    longitude of grid cell center

    units :

    degrees_east

    array([[189.22293223, 189.38990916, 189.55836619, ..., 293.77906088,
            294.0279241 , 294.27439931],
           [188.96836986, 189.13470591, 189.30253733, ..., 294.05584005,
            294.30444387, 294.55065969],
           [188.71234264, 188.87800731, 189.04515208, ..., 294.335053  ,
            294.58337453, 294.8292928 ],
           ...,
           [124.04724025, 123.88362026, 123.71852016, ...,  16.83171831,
             16.58436953,  16.33949649],
           [123.78686428, 123.62254238, 123.45672512, ...,  17.11814486,
             16.87043749,  16.62518298],
           [123.52798366, 123.36295986, 123.1964407 , ...,  17.40209947,
             17.1540526 ,  16.90845095]])

  • yc
    (y, x)
    float64
    16.53 16.78 17.02 ... 27.76 27.51

    long_name :

    latitude of grid cell center

    units :

    degrees_north

    array([[16.53498637, 16.7784556 , 17.02222429, ..., 27.36301592,
            27.11811045, 26.87289026],
           [16.69397341, 16.93865381, 17.18364512, ..., 27.5847719 ,
            27.33821848, 27.0913656 ],
           [16.85219179, 17.09808909, 17.34430872, ..., 27.80584314,
            27.55764558, 27.30915621],
           ...,
           [17.31179033, 17.56124674, 17.81104646, ..., 28.4502485 ,
            28.19718339, 27.94384744],
           [17.15589701, 17.40414034, 17.65272318, ..., 28.23129632,
            27.97989251, 27.72821596],
           [16.99919497, 17.24622904, 17.49358736, ..., 28.01160028,
            27.76185586, 27.51182726]])

• Data variables: (1)
  • Tair
    (time, y, x)
    float64
    nan nan nan ... 28.66 28.19 28.21

    units :

    C

    long_name :

    Surface air temperature

    type_preferred :

    double

    time_rep :

    instantaneous

    array([[[        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
            [        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
            [        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
            ...,
            [        nan,         nan,         nan, ..., 27.03290153,
             27.03125761, 27.33531541],
            [        nan,         nan,         nan, ..., 27.2784053 ,
             26.80261869, 27.08603517],
            [        nan,         nan,         nan, ..., 27.02344402,
             26.56473862, 26.73064933]],
    
           [[        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
            [        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
            [        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
    ...
            [        nan,         nan,         nan, ..., 27.8597472 ,
             27.82928439, 28.09249224],
            [        nan,         nan,         nan, ..., 27.89704094,
             27.31104941, 27.67387171],
            [        nan,         nan,         nan, ..., 27.46837113,
             27.0088944 , 27.23017976]],
    
           [[        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
            [        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
            [        nan,         nan,         nan, ...,         nan,
                     nan,         nan],
            ...,
            [        nan,         nan,         nan, ..., 28.95929072,
             28.87672039, 29.04890862],
            [        nan,         nan,         nan, ..., 29.036132  ,
             28.42273578, 28.68721201],
            [        nan,         nan,         nan, ..., 28.66381585,
             28.18595533, 28.20753022]]])

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(CFTimeIndex([1980-09-16 12:00:00, 1980-10-17 00:00:00, 1980-11-16 12:00:00,
                 1980-12-17 00:00:00, 1981-01-17 00:00:00, 1981-02-15 12:00:00,
                 1981-03-17 00:00:00, 1981-04-16 12:00:00, 1981-05-17 00:00:00,
                 1981-06-16 12:00:00, 1981-07-17 00:00:00, 1981-08-17 00:00:00,
                 1981-09-16 12:00:00, 1981-10-17 00:00:00, 1981-11-16 12:00:00,
                 1981-12-17 00:00:00, 1982-01-17 00:00:00, 1982-02-15 12:00:00,
                 1982-03-17 00:00:00, 1982-04-16 12:00:00, 1982-05-17 00:00:00,
                 1982-06-16 12:00:00, 1982-07-17 00:00:00, 1982-08-17 00:00:00,
                 1982-09-16 12:00:00, 1982-10-17 00:00:00, 1982-11-16 12:00:00,
                 1982-12-17 00:00:00, 1983-01-17 00:00:00, 1983-02-15 12:00:00,
                 1983-03-17 00:00:00, 1983-04-16 12:00:00, 1983-05-17 00:00:00,
                 1983-06-16 12:00:00, 1983-07-17 00:00:00, 1983-08-17 00:00:00],
                dtype='object', length=36, calendar='noleap', freq=None))

• Attributes: (11)

  title :

  /workspace/jhamman/processed/R1002RBRxaaa01a/lnd/temp/R1002RBRxaaa01a.vic.ha.1979-09-01.nc

  institution :

  U.W.

  source :

  RACM R1002RBRxaaa01a

  output_frequency :

  daily

  output_mode :

  averaged

  convention :

  CF-1.4

  references :

  Based on the initial model of Liang et al., 1994, JGR, 99, 14,415- 14,429.

  comment :

  Output from the Variable Infiltration Capacity (VIC) model.

  nco_openmp_thread_number :

  1

  NCO :

  netCDF Operators version 4.7.9 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)

  history :

  Fri Aug 7 17:57:38 2020: ncatted -a bounds,,d,, rasm.nc Tue Dec 27 14:15:22 2016: ncatted -a dimensions,,d,, rasm.nc rasm.nc Tue Dec 27 13:38:40 2016: ncks -3 rasm.nc rasm.nc history deleted for brevity

Just like the rectilinear grids GeoViews understands this kind of data natively. So we again wrap this dataset in a gv.Dataset and define a fixed range for the air teperature (Tair) values:

gvds = gv.Dataset(ds).redim.range(Tair=(-25, 25))
gvds

:Dataset   [time,xc,yc]   (Surface air temperature)

Now we can plot this data directly as a gv.QuadMesh, however this is generally quite slow, especially when we are working with bokeh where each grid point is rendered as a distinct polygon. We will therefore downsample the data by a factor of 3 along both dimensions:

opts.defaults(opts.Image(cmap='RdBu_r'), opts.QuadMesh(cmap='RdBu_r'))

quadmeshes = gvds.to(gv.QuadMesh, ['xc', 'yc'], dynamic=True)
quadmeshes.apply(lambda x: x.clone(x.data.Tair[::3, ::3])) * gv.feature.coastline

The problem is less severe when plotting using matplotlib but even then plotting can be fairly slow given a large enough mesh.

gv.output(quadmeshes.opts(cmap='RdBu_r') * gv.feature.coastline, backend='matplotlib', size=300)

If we want to explore a very large grid it therefore often makes sense to resample the data onto a rectilinear grid, which can be rendered much more efficiently. Once again we have the option of using the datashader based approach or the more accurate xESMF based approach.

Datashader

To regrid a QuadMesh using GeoViews we can import the rasterize operation. In the background the operation will convert the QuadMesh into a TriMesh, which datashader understands. To optimize this conversion so it occurs only when aggregating the QuadMesh for the first time we can activate the precompute option. Additionally we have to define an aggregator, in this case to compute the mean Tair value in a pixel:

from holoviews.operation.datashader import rasterize
rasterize(gv.project(quadmeshes), precompute=True, aggregator=dsh.mean('Tair')) * gv.feature.coastline

xESMF

Now we will once again use the xESMF based regridding for which we can still use the weighted_regrid operation, since it supports both rectilinear and curvilinear grids. Since the original data doesn’t have a very high resolution we will also disable the streams linking the operation to the plot dimensions and axis ranges.

from geoviews.operation.regrid import weighted_regrid
weighted_regrid(quadmeshes, streams=[]) * gv.feature.coastline

Finally lets clean up after ourselves one last time:

weighted_regrid.clean_weight_files()

Deleted 0 weight files

This web page was generated from a Jupyter notebook and not all interactivity will work on this website. Right click to download and run locally for full Python-backed interactivity.

---

Right click to download this notebook from GitHub.