{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Compare Mercator to ARGO\n\nThis notebook, we show:\n\n* how to easily rename variables and coordinates in an ARGO and a Mercator datasets,\n* how to plot a T-S diagram,\n* how to stack geographical coordinates,\n* how to compute a geographical distance between two datasets with coordinates,\n* how to interpolate a vertical profile with different methods.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Initialisations\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import xarray as xr\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport cartopy.crs as ccrs\nimport gsw\nimport cmocean\n\nimport xoa\nimport xoa.grid as xgrid\nimport xoa.regrid as xregrid\nimport xoa.cf as xcf\nimport xoa.coords as xcoords\nimport xoa.geo as xgeo\nimport xoa.plot as xplot\n\nmpl.rc(\"axes\", grid=True)\nxr.set_options(display_style=\"text\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Register the `xoa <accessors>` accessors :\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xoa.register_accessors()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Register the Mercator and ARGO naming specifications\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xcf.register_cf_specs(xoa.get_data_sample(\"mercator.cfg\"))\nxcf.register_cf_specs(xoa.get_data_sample(\"argo.cfg\"))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here is what these CF specifications contain\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "with open(xoa.get_data_sample(\"mercator.cfg\")) as f:\n    print(f.read())\nwith open(xoa.get_data_sample(\"argo.cfg\")) as f:\n    print(f.read())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read and decode the datasets\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### ARGO profiles\n\nThe ARGO data were downloaded with the `argopy <https://argopy.readthedocs.io/en/latest/>`_\npackage with a piece of code close to the following::\n\n    from argopy import DataFetcher\n    adf = DataFetcher().float(7900573)\n    ds = adf.to_xarray().argo.point2profile().isel(N_PROF=slice(-10, None))\n    ds.to_netcdf(\"argo-7900573.nc\")\n\nHere is the web page of this float: https://fleetmonitoring.euro-argo.eu/float/7900573\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ds_argo = xoa.open_data_sample(\"argo-7900573.nc\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We decode the names to make them more generic thanks to the CF specifications\nand select only the variables of interest.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ds_argo = ds_argo.xoa.decode()[[\"temp\", \"sal\", \"pres\"]]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We compute depths with the :mod:`gsw` package and assign them\nas coordinates of the ARGO dataset.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "alat = ds_argo.lat.broadcast_like(ds_argo.pres)\nadepth = gsw.z_from_p(ds_argo.pres, alat)\nadepth = adepth.where(adepth.notnull(), adepth.min())\nadepth.attrs.update({\"long_name\": \"Depth\", \"units\": \"m\"})\nds_argo = ds_argo.assign_coords(depth=adepth)\nds_argo = ds_argo.isel(level=slice(None, None, -1))\nds_argo = ds_argo.set_index(N_PROF=\"time\").rename(N_PROF=\"time\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Quick plot of the salinity that highligths the mediterranean water.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(ncols=2, figsize=(14, 6))\nds_argo.sal.plot.contourf(\"time\", \"depth\", cmap=\"cmo.haline\", levels=20, ax=axs[0])\nds_argo.sal.plot.contour(\"time\", \"depth\", levels=[35.65], linewidths=1, colors=\"k\", ax=axs[0])\nxplot.plot_ts(\n    ds_argo.temp,\n    ds_argo.sal,\n    potential=False,\n    scatter_c=ds_argo.depth,\n    scatter_cmap=\"cmo.deep\",\n    axes=axs[1],\n)\naxs[1].axvline(x=35.65, ls=\"--\", color=\"k\");"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we select the last ARGO profile as a reference profile\nsince it better samples de mediterranean water vein.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ds_argo_prof = ds_argo.isel(time=-1).squeeze(drop=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Mercator\n\nThe Mercator data come from the IBI configuration and were downloaded\nfrom the CMEMS site:\nhttps://resources.marine.copernicus.eu/product-detail/IBI_ANALYSISFORECAST_PHY_005_001/INFORMATION\nThe dataset has been undersampled by a factor 3 in longitude and latitude\nto limit disk usage.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ds_merc = xoa.open_data_sample(\"ibi-argo-7900573.nc\").xoa.decode()\nds_merc = ds_merc.isel(depth=slice(None, None, -1))\nds_merc = ds_merc.assign_coords(depth=-ds_merc.depth)\nprint(ds_merc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "By analogy with the ensemble data assimilation,\nthis block of data can be viewed as a ensemble of profiles that emulates\nthe uncertainties of the model, so we stack geographical coordinates in a `member` dimension\nthanks to :func:`~xoa.coords.geo_stack`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ds_merc_ens_rect = xcoords.geo_stack(ds_merc, \"member\")\nprint(ds_merc_ens_rect)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is a rectangular selection and we find it more appropriate to retain\nonly profiles that are in a given radius of proximity, that is chosen\nclose to the surface radius of deformation.\nWe compute the distance from each model profile to the selected ARGO profile\nwith :func:`~xoa.geo.get_distances`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "radius = 35e3  # 35 km\ndist_merc2argo = (\n    xgeo.get_distances(ds_merc_ens_rect, ds_argo_prof)\n    .squeeze()\n    .rename(npts0=\"member\")\n    .assign_coords(member=ds_merc_ens_rect.member)\n)\nds_merc_ens = ds_merc_ens_rect.where(dist_merc2argo < radius, drop=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compare Mercator and ARGO\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot the situation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pmerc = ccrs.Mercator()\npcar = ccrs.PlateCarree()\nfig, ax = plt.subplots(subplot_kw={\"projection\": pmerc}, figsize=(8, 8))\nax.set_extent(xgeo.get_extent(ds_merc, margin=1, square=True))\nax.gridlines(draw_labels=True)\nax.add_wms(\"https://ows.emodnet-bathymetry.eu/wms\", \"emodnet:mean_atlas_land\", alpha=0.5)\nkws = dict(c=\"C3\", s=15, marker=\"s\", transform=pcar)\nax.scatter(ds_merc_ens_rect.lon, ds_merc_ens_rect.lat, label=\"Rejected\", alpha=0.15, **kws)\nax.scatter(ds_merc_ens.lon, ds_merc_ens.lat, label='Mercator', **kws)\nax.scatter(ds_argo_prof.lon, ds_argo_prof.lat, s=100, c=\"C2\", transform=pcar, label='ARGO')\nplt.legend();"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We used the :func:`xoa.geo.get_extent` to compute a square geographical extent based\non the coordinates of a dataset with an added margin.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now interpolate Mercator onto ARGO time and position for both the ensemble and\nthe reference profile. All interpolations can be performed\nwith the :meth:`~xarray.Dataset.interp` method since we are dealing with\naxis (1D) coordinates.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ds_merc_ens = ds_merc_ens.interp(time=ds_argo_prof.time.values)\nds_merc_ens_std = ds_merc_ens.std(\"member\")\nds_merc_prof = ds_merc.interp(time=ds_argo_prof.time, lon=ds_argo_prof.lon, lat=ds_argo_prof.lat)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We plot now the model, its simulated uncertainty and the observations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(4, 7))\nplt.fill_betweenx(\n    ds_merc_prof.sal.depth,\n    ds_merc_prof.sal - 1.96 * ds_merc_ens_std.sal,\n    ds_merc_prof.sal + 1.96 * ds_merc_ens_std.sal,\n    fc=\"C1\",\n    alpha=0.2,\n)\nds_merc_prof.sal.plot(y=\"depth\", color=\"C1\", label=\"Mercator\")\nplt.plot(ds_merc_prof.sal.values, ds_merc_prof.depth, \"o-\", color=\"C1\")\nds_argo_prof.sal.plot(y=\"depth\", label=\"ARGO\", color=\"C0\")\nplt.legend();"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Therefore, if we accept an uncertainty in the positioning of the ocean\nstructure in the model, it differences with (assimilated) observations\nbecome acceptables.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As you can see on the previous figures, the ARGO profiles come with\na very high vertical resolution at all depths, whereas the model\ncomes with a crude vertical resolution far from the surface.\nEach model thermodynamical model value is representative of the\ncube that defines its XYZ cell.\nThis means that despite the model does simulate the small scale variability,\nthe cell average may be accurate.\nSo we propose here to perform comparisons in the model\nand the observational spaces. We use the function :func:`xoa.regrid.regrid1d`\nthat supports the linear (:attr:`~xoa.regrid.regrid1d_methods.linear` method) interpolation\nto higher resolutions and the cell averaging\n(:attr:`~xoa.regrid.regrid1d_methods.cellave` method) to lower resolutions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sal_merc_ens_o = xregrid.regrid1d(ds_merc_ens.sal, ds_argo_prof.depth, method=\"linear\")\nsal_merc_prof_o = xregrid.regrid1d(ds_merc_prof.sal, ds_argo_prof.depth, method=\"linear\")\nsal_argo_prof_m = xregrid.regrid1d(ds_argo_prof.sal, ds_merc_prof.depth, method=\"cellave\")\nmerc_prof_edges = xgrid.get_edges(ds_merc_prof.depth, \"depth\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot them.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=(10, 5))\nsal_merc_prof_o.plot(y=\"depth\", color=\"C1\", label=\"Mercator\", ax=axes[0])\nds_argo_prof.sal.plot(y=\"depth\", label=\"ARGO\", color=\"C0\", ax=axes[0])\naxes[0].set_title(\"On ARGO depths\")\naxis = axes[0].axis()\nds_merc_prof.sal.plot(y=\"depth\", color=\"C1\", label=\"Mercator\", ax=axes[1])\nplt.stairs(ds_merc_prof.sal, merc_prof_edges, orientation=\"horizontal\", color=\"C1\", lw=0.5)\nsal_argo_prof_m.plot(y=\"depth\", label=\"ARGO\", color=\"C0\", ax=axes[1])\nplt.stairs(sal_argo_prof_m, merc_prof_edges, orientation=\"horizontal\", color=\"C0\", lw=0.5)\naxes[1].set_title(\"On Mercator depths\")\naxes[1].axis(axis);"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "These plots show that the model is too haline, except at the\nmediterranean water depth where the model is very close to the observations\nonce an appropriate vertical regridding is performed.\nHowever, the vertical resolution is not sufficient enough to simulate\nthe maximum of salinity.\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.10"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}