xbatcher for Machine Learning Part 1


Here, we will be covering how to use xbatcher with Keras/Tensorflow convolutional neural network (CNN) models.








Strongly Recommended

Not strictly needed to understand this tutorial

This notebook replicates the work of Sinha and Abernathey, 2021, where the goal is to use a CNN to learn ocean surface currents (which are usually inferred diagnostically or modelled) from variables that can be observed directly, like sea surface temperature (SST) or wind stress.

Can we learn to predict ocean currents with just one snapshot of data?


To start, let’s import some libraries we’ll need. The important libraries here are numpy, xarray, xbatcher and tensorflow, while most of the others aren’t strictly necessary.

import numpy as np
import xarray as xr

from intake import open_catalog
from dataclasses import dataclass
from typing import Iterable
from matplotlib import pyplot as plt
from IPython.display import clear_output
import tensorflow as tf
2023-12-08 14:40:52.461453: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: SSE4.1 SSE4.2 AVX AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
import xbatcher as xb

Designing Scenarios

We want to experiment with different neural network models by providing different inputs, and perhaps by playing with whether or not we run them through a convolutional layer. There are a lot of possibilities here, and if we approach it haphazardly, we’ll end up with a mess of scattered experiments and results mixed in with other code.

Instead, we can be more systematic about it. We know we want to define an individual scenario once, and then have it stay constant through the workflow. This way, there will be no complexities later on about whether we’re referring to the right dataset, etc. With that in mind, we should use a dataclass. We want something minimal here, just enough to store the names of variables we’re interested in.

What is the structure of each experiment? We want some input variables to be run through a 2D convolutional layer, while some other inputs will be passed through directly to the dense part of the neural network. Both of these can be lists of strings, so we define conv_var and input_var as Iterable[str].

Likewise, we have more than one target, so we define the target item as Iterable[str] as well. Outside of the Scenario dataclass, we define target as a list: ['U', 'V']. Since we’re only interested in learning the currents, this won’t change.

Finally, we need to name each scenario something distinct, so when we create data subsets for training, testing, and prediction, we can recover them later.

class Scenario:
    conv_var: Iterable[str]
    input_var: Iterable[str]
    target: Iterable[str]
    name: str
target = ['U', 'V']
sc1 = Scenario(['SSH'],             ['TAUX', 'TAUY'], target, name = "derp")
sc5 = Scenario(['SSH', 'SST'], ['X', 'TAUX', 'TAUY'], target, name = "herp")

Data and Preprocessing

For our dataset, we will be using ocean data from a high-resolution CESM POP model.

We have some necessary I/O routines, but they aren’t central to our problem, aside from the addtion of the new variables X, Y, Z, dx and dy, which represent Euclidean positions and distances between grid points.

You can have a look in the notebook below if you’re curious about it.

%run ./surface_currents_prep.ipynb

From this notebook, we get a few new functions.

  • prepare_data takes a scenario, as well as the time slices for training, testing, and prediction we are interested in, and the time slice we’ll use for the NaN mask. It adds the new grid variables, and then stores each slice in a new zarr store that we can access later. This speeds up future I/O, which is helpful when modifying the model. Each scenario is stored separately.

  • load_training_data loads the training data created for the scenario passed to it.

  • load_test_data loads the testing data created for the scenario passed to it.

  • load_predict_data loads the prediction input data created for the scenario passed to it.

You can comment out prepare_data after you’ve run it once, it will save time if you rerun the whole notebook again.

prepare_data(sc5, 200, 1000, 1000, 200)
ValueError                                Traceback (most recent call last)
Cell In[9], line 1
----> 1 prepare_data(sc5, 200, 1000, 1000, 200)

File /tmp/ipykernel_3348/3337232756.py:3, in prepare_data(sc, training_time, test_time, predict_time, mask_time)
      1 def prepare_data(sc, training_time, test_time, predict_time, mask_time=11):
      2     cat = open_catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/ocean/CESM_POP.yaml")
----> 3     ds  = cat["CESM_POP_hires_control"].to_dask()
      4     ds = ds.rename({'U1_1':'U', 'V1_1':'V', 'TAUX_2':'TAUX', 'TAUY_2':'TAUY', 'SSH_2':'SSH', 'ULONG':'XU', 'ULAT':'YU'})
      5     ds = add_grid(ds)

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/intake_xarray/base.py:69, in DataSourceMixin.to_dask(self)
     67 def to_dask(self):
     68     """Return xarray object where variables are dask arrays"""
---> 69     return self.read_chunked()

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/intake_xarray/base.py:44, in DataSourceMixin.read_chunked(self)
     42 def read_chunked(self):
     43     """Return xarray object (which will have chunks)"""
---> 44     self._load_metadata()
     45     return self._ds

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/intake/source/base.py:283, in DataSourceBase._load_metadata(self)
    281 """load metadata only if needed"""
    282 if self._schema is None:
--> 283     self._schema = self._get_schema()
    284     self.dtype = self._schema.dtype
    285     self.shape = self._schema.shape

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/intake_xarray/base.py:18, in DataSourceMixin._get_schema(self)
     15 self.urlpath = self._get_cache(self.urlpath)[0]
     17 if self._ds is None:
---> 18     self._open_dataset()
     20     metadata = {
     21         'dims': dict(self._ds.dims),
     22         'data_vars': {k: list(self._ds[k].coords)
     23                       for k in self._ds.data_vars.keys()},
     24         'coords': tuple(self._ds.coords.keys()),
     25     }
     26     if getattr(self, 'on_server', False):

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/intake_xarray/xzarr.py:46, in ZarrSource._open_dataset(self)
     44     self._ds = xr.open_mfdataset(self.urlpath, **kw)
     45 else:
---> 46     self._ds = xr.open_dataset(self.urlpath, **kw)

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/xarray/backends/api.py:572, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, inline_array, chunked_array_type, from_array_kwargs, backend_kwargs, **kwargs)
    560 decoders = _resolve_decoders_kwargs(
    561     decode_cf,
    562     open_backend_dataset_parameters=backend.open_dataset_parameters,
    568     decode_coords=decode_coords,
    569 )
    571 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 572 backend_ds = backend.open_dataset(
    573     filename_or_obj,
    574     drop_variables=drop_variables,
    575     **decoders,
    576     **kwargs,
    577 )
    578 ds = _dataset_from_backend_dataset(
    579     backend_ds,
    580     filename_or_obj,
    590     **kwargs,
    591 )
    592 return ds

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/xarray/backends/zarr.py:992, in ZarrBackendEntrypoint.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, synchronizer, consolidated, chunk_store, storage_options, stacklevel, zarr_version)
    971 def open_dataset(  # type: ignore[override]  # allow LSP violation, not supporting **kwargs
    972     self,
    973     filename_or_obj: str | os.PathLike[Any] | BufferedIOBase | AbstractDataStore,
    989     zarr_version=None,
    990 ) -> Dataset:
    991     filename_or_obj = _normalize_path(filename_or_obj)
--> 992     store = ZarrStore.open_group(
    993         filename_or_obj,
    994         group=group,
    995         mode=mode,
    996         synchronizer=synchronizer,
    997         consolidated=consolidated,
    998         consolidate_on_close=False,
    999         chunk_store=chunk_store,
   1000         storage_options=storage_options,
   1001         stacklevel=stacklevel + 1,
   1002         zarr_version=zarr_version,
   1003     )
   1005     store_entrypoint = StoreBackendEntrypoint()
   1006     with close_on_error(store):

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/xarray/backends/zarr.py:464, in ZarrStore.open_group(cls, store, mode, synchronizer, group, consolidated, consolidate_on_close, chunk_store, storage_options, append_dim, write_region, safe_chunks, stacklevel, zarr_version, write_empty)
    461             raise FileNotFoundError(f"No such file or directory: '{store}'")
    462 elif consolidated:
    463     # TODO: an option to pass the metadata_key keyword
--> 464     zarr_group = zarr.open_consolidated(store, **open_kwargs)
    465 else:
    466     zarr_group = zarr.open_group(store, **open_kwargs)

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/zarr/convenience.py:1351, in open_consolidated(store, metadata_key, mode, **kwargs)
   1348         metadata_key = "meta/root/consolidated/" + metadata_key
   1350 # setup metadata store
-> 1351 meta_store = ConsolidatedStoreClass(store, metadata_key=metadata_key)
   1353 # pass through
   1354 chunk_store = kwargs.pop("chunk_store", None) or store

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/zarr/storage.py:2947, in ConsolidatedMetadataStore.__init__(self, store, metadata_key)
   2944 self.store = Store._ensure_store(store)
   2946 # retrieve consolidated metadata
-> 2947 meta = json_loads(self.store[metadata_key])
   2949 # check format of consolidated metadata
   2950 consolidated_format = meta.get("zarr_consolidated_format", None)

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/zarr/storage.py:1428, in FSStore.__getitem__(self, key)
   1426 key = self._normalize_key(key)
   1427 try:
-> 1428     return self.map[key]
   1429 except self.exceptions as e:
   1430     raise KeyError(key) from e

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/fsspec/mapping.py:151, in FSMap.__getitem__(self, key, default)
    149 k = self._key_to_str(key)
    150 try:
--> 151     result = self.fs.cat(k)
    152 except self.missing_exceptions:
    153     if default is not None:

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/fsspec/asyn.py:118, in sync_wrapper.<locals>.wrapper(*args, **kwargs)
    115 @functools.wraps(func)
    116 def wrapper(*args, **kwargs):
    117     self = obj or args[0]
--> 118     return sync(self.loop, func, *args, **kwargs)

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/fsspec/asyn.py:103, in sync(loop, func, timeout, *args, **kwargs)
    101     raise FSTimeoutError from return_result
    102 elif isinstance(return_result, BaseException):
--> 103     raise return_result
    104 else:
    105     return return_result

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/fsspec/asyn.py:56, in _runner(event, coro, result, timeout)
     54     coro = asyncio.wait_for(coro, timeout=timeout)
     55 try:
---> 56     result[0] = await coro
     57 except Exception as ex:
     58     result[0] = ex

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/fsspec/asyn.py:446, in AsyncFileSystem._cat(self, path, recursive, on_error, batch_size, **kwargs)
    444     ex = next(filter(is_exception, out), False)
    445     if ex:
--> 446         raise ex
    447 if (
    448     len(paths) > 1
    449     or isinstance(path, list)
    450     or paths[0] != self._strip_protocol(path)
    451 ):
    452     return {
    453         k: v
    454         for k, v in zip(paths, out)
    455         if on_error != "omit" or not is_exception(v)
    456     }

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/asyncio/tasks.py:408, in wait_for(fut, timeout)
    405 loop = events.get_running_loop()
    407 if timeout is None:
--> 408     return await fut
    410 if timeout <= 0:
    411     fut = ensure_future(fut, loop=loop)

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/gcsfs/core.py:1027, in GCSFileSystem._cat_file(self, path, start, end, **kwargs)
   1025 else:
   1026     head = {}
-> 1027 headers, out = await self._call("GET", u2, headers=head)
   1028 return out

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/gcsfs/core.py:437, in GCSFileSystem._call(self, method, path, json_out, info_out, *args, **kwargs)
    433 async def _call(
    434     self, method, path, *args, json_out=False, info_out=False, **kwargs
    435 ):
    436     logger.debug(f"{method.upper()}: {path}, {args}, {kwargs.get('headers')}")
--> 437     status, headers, info, contents = await self._request(
    438         method, path, *args, **kwargs
    439     )
    440     if json_out:
    441         return json.loads(contents)

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/decorator.py:221, in decorate.<locals>.fun(*args, **kw)
    219 if not kwsyntax:
    220     args, kw = fix(args, kw, sig)
--> 221 return await caller(func, *(extras + args), **kw)

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/gcsfs/retry.py:123, in retry_request(func, retries, *args, **kwargs)
    121     if retry > 0:
    122         await asyncio.sleep(min(random.random() + 2 ** (retry - 1), 32))
--> 123     return await func(*args, **kwargs)
    124 except (
    125     HttpError,
    126     requests.exceptions.RequestException,
    129     aiohttp.client_exceptions.ClientError,
    130 ) as e:
    131     if (
    132         isinstance(e, HttpError)
    133         and e.code == 400
    134         and "requester pays" in e.message
    135     ):

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/gcsfs/core.py:430, in GCSFileSystem._request(self, method, path, headers, json, data, *args, **kwargs)
    427 info = r.request_info  # for debug only
    428 contents = await r.read()
--> 430 validate_response(status, contents, path, args)
    431 return status, headers, info, contents

File ~/miniconda3/envs/xbatcher-ML-1-cookbook-dev/lib/python3.10/site-packages/gcsfs/retry.py:108, in validate_response(status, content, path, args)
    106     raise requests.exceptions.ProxyError()
    107 elif "invalid" in str(msg):
--> 108     raise ValueError(f"Bad Request: {path}\n{msg}")
    109 elif error:
    110     raise HttpError(error)

ValueError: Bad Request: https://storage.googleapis.com/download/storage/v1/b/pangeo-cesm-pop/o/control%2F.zmetadata?alt=media
User project specified in the request is invalid.

Next, we’ll load our training data and pick out the part we want to train with.

NOTE: Coordinates and attributes are dropped for speed, doing this shouldn’t be necessary in future (optimized) versions of xarray/xbatcher.

ds_training = load_training_data(sc5)
ds_training = just_the_data(ds_training)

Looking inside ds_training, we see only the variables we would expect from sc5.

ds_training = select_from(ds_training)

Model Setup

We have a model architecuture we’re happy with already defined, so for this tutorial, we’ll focus on how to use xbatcher to generate training sets for the model. From the notebook below, we recieve:

  • get_model() Creates a mixed neural network based on some parameters. The architecture is intentionally a little arbitrary in terms of the depth of the dense part of the network, the depth of the convolutional part of the network, and the convolution kernel size. Returns a compiled Keras model.

  • LossHistory() Only needed here because it has to be passed to model.fit().

  • train() We will walk through this routine below.

Have a look inside for more details!

%run ./surface_currents_model.ipynb

Now the fun part: we define the train function to deal with high-level aspects of training the model, which means this is a good place to use xbatcher. Let’s walk through it…

The arguments to train are

  • ds: xr.DataSet The dataset you want to work with.

  • sc: Scenario The scenario you want to work with.

  • conv_dims: List[int] This is the shape of the stencil that will be passed to the first convolutional layer. We are only interested in 2D convolutions here, so it will need to be a list of two integers. Note that this is distinct from the convolutional kernel.

  • nfilters: int How many filters do we want to map the first convolution layer to?

  • conv_kernels: List[int] Each entry denotes the convolution kernel of a new convolution layer. train works best for odd-numbered convolution kernels.

  • dense_layers: int The number of dense layers in the model.

For this example, we only use one convolution layer, which makes some things simpler. Feel free to experiment with these parameters to use different data sets and create new CNN models.

sc = sc5
conv_dims = [5,5]
nfilters = 80
conv_kernels = [5]
dense_layers = 3

We’ll need some info about how to rectify the output of the convolution layers with raw input from other variables (see the surface_currents_model.ipynb notebook for more info). Based on the convolution kernel, we know how the output of a convolution layer will be shaped compared to the input: a halo of a certain size will be removed from the edges. For odd convolution kernels, the halo thickness is always \(\frac{n - 1}{2}\) where \(n\) is the kernel.

halo_size = int((np.sum(conv_kernels) - len(conv_kernels))/2)

Training a Model with xbatcher

Since we are trying to learn from a single 2D snapshot, it makes sense to iterate in both latitude and longitude. What we want are individual samples of the size given by conv_dims, but batched in a way that we can pass the correct number of samples to the model as a single tensor. So, input_dims will contain entries for both nlon and nlat. To take full advantage of the available data, we can add an overlap to make sure halo points are fully included in the neighboring samples.

NOTE: xbatcher currently runs slowly with concat_input_dims=True, and running without it will result in batches of size one. Therefore, we use an implemenatation of xarray rolling to mimic what xbatcher does. This is not good strategy when using large datasets, but for this example, the differences are minimal. We anticipate that fixed-size batches and some optimizations will be implemented in xbatcher in the future.

nlons, nlats = conv_dims
# bgen = xb.BatchGenerator(
#     ds_training,
#     {'nlon':nlons,       'nlat':nlats},
#     {'nlon':2*halo_size, 'nlat':2*halo_size}
# )
latlen = len(ds_training['nlat'])
lonlen = len(ds_training['nlon'])
nlon_range = range(nlons,lonlen,nlons - 2*halo_size)
nlat_range = range(nlats,latlen,nlats - 2*halo_size)

batch = (
    .rolling({"nlat": nlats, "nlon": nlons})
    .construct({"nlat": "nlat_input", "nlon": "nlon_input"})[{'nlat':nlat_range, 'nlon':nlon_range}]
    .stack({"input_batch": ("nlat", "nlon")}, create_index=False)
    .rename_dims({'nlat_input':'nlat', 'nlon_input':'nlon'})
    # .chunk({'input_batch':32, 'nlat':nlats, 'nlon':nlons})
rnds = list(range(len(batch['input_batch'])))
batch = batch[{'input_batch':(rnds)}]
# use with rolling
def batch_generator(batch_set, batch_size):
    n = 0
    while n < len(batch_set['input_batch']) - batch_size:
        yield batch_set.isel({'input_batch':range(n,(n+batch_size))})
        n += batch_size
# # use with xbatcher
# def batch_generator(bgen, batch_size):
#     b = (batch for batch in bgen)
#     n = 0
#     while n < 400:
#         batch_stack = [ next(b) for i in range(batch_size) ]
#         yield xr.concat(batch_stack, 'sample')
#         n += 1
bgen = batch_generator(batch, 4096)
# bgen = batch_generator(bgen, 32)

We need a subsetting stencil (sub) to compensate for the fact that a halo is removed by each convolution layer. This means that the input_var variables will be the wrong size at the concat layer unless we strip the halo from them.

sub = {'nlon':range(halo_size,nlons-halo_size),

Here, we generate our model and our history callback.

model = get_model(halo_size, ds_training, sc, conv_dims, nfilters, conv_kernels, dense_layers)
history = LossHistory()

And now, we can construct our training loop. Most use cases of the xb.BatchGenerator will take the form of a for-loop with the construct for batch in bgen.

Once we have a batch, we still have some things to do before we can pass the data to the model.

So when we look at the contents of each batch, we see

# a = []
# for batch in bgen:
#     a = batch
#     break
# a

…but our model expects tensors where the different variables are stacked in a new dimension we will call var.

Looking at model.fit(), we have two separate inputs because of the distinction between convolved inputs and raw inputs. Therefore, the model expects these inputs to be given as a list of the two. The training target is relatively straightforward. On the next line, we have a couple of parameters we can experiment with. The important thing to note is the batch_size parameter; you may need to check that the sample dimension is compatible with the dimensions that xb.BatchGenerator returned. And finally, we pass our history class as a callback so we can see how the model training is progressing.

for batch in bgen:
    batch_conv   = [batch[x] for x in sc.conv_var]
    batch_input  = [batch[x][sub] for x in sc.input_var]
    batch_target = [batch[x][sub] for x in sc.target]
    batch_conv   = xr.merge(batch_conv).to_array('var').transpose(...,'var')
    batch_input  = xr.merge(batch_input).to_array('var').transpose(...,'var')
    batch_target = xr.merge(batch_target).to_array('var').transpose(...,'var')

    model.fit([batch_conv, batch_input],
              batch_size=32, verbose=0,# epochs=4,

And now that we have our model trained, we can save it for future use. Note that once this model is saved, we don’t need to rerun much from above to continue with testing or prediction.

model.save('models/'+ sc.name)
np.savez('models/history_'+sc.name, losses=history.mae, mse=history.mse, accuracy=history.accuracy)

Training Function

#train(ds_training, sc5, conv_dims, conv_kernels)

Testing the Model

ds_test = load_test_data(sc5)
ds_test = just_the_data(ds_test)
ds_test = select_from(ds_test)
latlen = len(ds_test['nlat'])
lonlen = len(ds_test['nlon'])
nlon_range = range(nlons,lonlen,nlons - 2*halo_size)
nlat_range = range(nlats,latlen,nlats - 2*halo_size)

batch_test = (
    .rolling({"nlat": nlats, "nlon": nlons})
    .construct({"nlat": "nlat_input", "nlon": "nlon_input"})[{'nlat':nlat_range, 'nlon':nlon_range}]
    .stack({"input_batch": ("nlat", "nlon")}, create_index=False)
    .rename_dims({'nlat_input':'nlat', 'nlon_input':'nlon'})
    # .chunk({'input_batch':32, 'nlat':nlats, 'nlon':nlons})

Let’s load the trained model from before:

model = tf.keras.models.load_model('models/'+ sc.name, custom_objects={'Grid_MAE':Grid_MAE})
test_conv   = [batch_test[x]      for x in sc.conv_var]
test_input  = [batch_test[x][sub] for x in sc.input_var]
test_target = [batch_test[x][sub] for x in sc.target]
test_conv   = xr.merge(test_conv  ).to_array('var').transpose(...,'var')
test_input  = xr.merge(test_input ).to_array('var').transpose(...,'var')
test_target = xr.merge(test_target).to_array('var').transpose(...,'var')
model.evaluate([test_conv, test_input], test_target)

Testing Function

#test(ds_test, sc5, conv_dims, conv_kernels)

Making Predictions

ds_predict = load_predict_data(sc5)
ds_predict = just_the_data(ds_predict)
ds_predict = select_from(ds_predict)
latlen = len(ds_predict['nlat'])
lonlen = len(ds_predict['nlon'])
nlon_range = range(nlons,lonlen,nlons - 2*halo_size)
nlat_range = range(nlats,latlen,nlats - 2*halo_size)

batch_predict = (
    .rolling({"nlat": nlats, "nlon": nlons})
    .construct({"nlat": "nlat_input", "nlon": "nlon_input"})[{'nlat':nlat_range, 'nlon':nlon_range}]
    .stack({"input_batch": ("nlat", "nlon")}, create_index=False)
    .rename_dims({'nlat_input':'nlat', 'nlon_input':'nlon'})
    # .chunk({'input_batch':32, 'nlat':nlats, 'nlon':nlons})
model = tf.keras.models.load_model('models/'+ sc.name, custom_objects={'Grid_MAE':Grid_MAE})
predict_conv  = [batch_test[x]      for x in sc.conv_var]
predict_input = [batch_test[x][sub] for x in sc.input_var]
predict_conv  = xr.merge(predict_conv ).to_array('var').transpose(...,'var')
predict_input = xr.merge(predict_input).to_array('var').transpose(...,'var')
predict_target = model.predict([predict_conv, predict_input])

Prediction Function

#predict_target = predict(ds_predict, sc5, conv_dims, conv_kernels)

Prediction Results

Now, let’s take a look at the predicted surface currents and see how the model did. Notice that the predicted data can be retrieved fairly easily from our default setup. We only have to reshape them, with respect to the original dimensions and a halo that will be stripped off. This is because we chose to make the convolution kernel equal to the dimensions of the samples, which means the model will give results at individual points.

However, the convolution kernal can be different, it’s just that we will then have to use a more complex process to restructure our grid.

Note also that if there were nans removed, we would have to keep track of how to map the unstructured model inputs back to the original grid and insert nans in the correct positions.

U = ds_predict['U']
V = ds_predict['V']
U_pred = predict_target[:,0,0,0].reshape(545, 345)
V_pred = predict_target[:,0,0,1].reshape(545, 345)
plt.figure(figsize=(10, 6))
plt.pcolormesh(U, cmap='RdBu_r')
plt.clim([-100, 100])
plt.figure(figsize=(10, 6))
plt.pcolormesh(U_pred, cmap='RdBu_r')
plt.clim([-100, 100])

We can see that they look very similar, but to get a better idea of what our errors look like, we can subtract them.

plt.figure(figsize=(10, 6))
plt.pcolormesh(U_pred - U[3:-2,3:-2], cmap='RdBu_r') # double-check U indexing
plt.clim([-100, 100])