How to use the Mercator OpenDap

The access to this opendap requires authentication. To register please contact the Mercator-ocean service desk by :

Email: products@mercator-ocean.fr

Mandatory pre-requisite

The LDAP-backed HTTP/S-Basic authentication should work by reading credentials from the .netrc file given that the .dodsrc file is set to point to them. Here's a short summary of the configuration Add your URS/LDAP credentials to the .netrc file, associating them with the URS/OpenDAP server that you normally authenticate with, like this:

machine tds.mercator-ocean.fr
                login <your login>
                password <your password>

Next, edit the .dodsrc file in your HOME directory so that it tells DAP clients to use the .netrc file for password information, the line have to look like :

HTTP.NETRC=/home/jdoe/.netrc

At the end, you must change the rights on the .netrc file using the following command:

chmod 600 .netrc

First part : Get data with python

To be able to use theses python scripts, you will need:

  • The Python3.x under anaconda
  • The following modules: numpy, pandas, xarray, pydap and requests
  • If you encounter timeout problems, please try to split your data retrieval requests over time

To install a module you can use the following command : conda install module_name

1 - Connect opendap

You have to select a dataset and get the opendap adress which look like : '' 'http://tds.mercator-ocean.fr/thredds/dodsC/psy4v3r1-hourly-gridT'''

Dataset described variables on an irregular grid : ORCA and Arakawa type C grid.
So we have to convert our longitude and latitude in indexes of array.

In [1]:
import numpy as np
import pandas as pd
import xarray as xr
In [2]:
datapgn = xr.open_dataset('https://tds.mercator-ocean.fr/thredds/dodsC/psy4v3r1-daily-gridT')
In [3]:
print(datapgn)
<xarray.Dataset>
Dimensions:       (deptht: 50, time_counter: 4844, x: 4322, y: 3059)
Coordinates:
    nav_lon       (y, x) float32 ...
    nav_lat       (y, x) float32 ...
  * deptht        (deptht) float32 0.49402538 1.5413754 ... 5274.784 5727.9165
  * x             (x) int32 1 2 3 4 5 6 7 ... 4316 4317 4318 4319 4320 4321 4322
  * y             (y) int32 1 2 3 4 5 6 7 ... 3053 3054 3055 3056 3057 3058 3059
  * time_counter  (time_counter) datetime64[ns] 2006-10-11T12:00:00 ... 2020-01-14T12:00:00
Data variables:
    votemper      (time_counter, deptht, y, x) float32 ...
Attributes:
    Conventions:  CF-1.0
    file_name:    ORCA12_LIM-T00_y2020m01d13_gridT.nc
    institution:  MERCATOR OCEAN
    source:       NEMO
    TimeStamp:    2020-JAN-19 20:05:20 GMT-0000
    references:   http://www.mercator-ocean.eu

2 - Get dataset information : dimensions and variables

In [3]:
datapgn
Out[3]:
<xarray.Dataset>
Dimensions:       (deptht: 50, time_counter: 4844, x: 4322, y: 3059)
Coordinates:
    nav_lon       (y, x) float32 ...
    nav_lat       (y, x) float32 ...
  * deptht        (deptht) float32 0.49402538 1.5413754 ... 5274.784 5727.9165
  * x             (x) int32 1 2 3 4 5 6 7 ... 4316 4317 4318 4319 4320 4321 4322
  * y             (y) int32 1 2 3 4 5 6 7 ... 3053 3054 3055 3056 3057 3058 3059
  * time_counter  (time_counter) datetime64[ns] 2006-10-11T12:00:00 ... 2020-01-14T12:00:00
Data variables:
    votemper      (time_counter, deptht, y, x) float32 ...
Attributes:
    Conventions:  CF-1.0
    file_name:    ORCA12_LIM-T00_y2020m01d13_gridT.nc
    institution:  MERCATOR OCEAN
    source:       NEMO
    TimeStamp:    2020-JAN-19 20:05:20 GMT-0000
    references:   http://www.mercator-ocean.eu

Area selection

We define area extension and then translate it in indexes of X and Y dimensions.
After we will use the method isel to extract data over our desired area.
The irregular grid is define with a first index "X" corresponding to a longitude equal to 73°. Moreover the two first value along x-axis are duplicate at the end of the array.
If you want to extract an area with longitude crossing 73° and not global, there is some tricks presented in part 4.

Time selection

For time, you can use index method (isel) or value method (sel).
You can combine method sel with argument nearest which give you the closest value to your request.
For multiple selection, you can use the method date_range from pandas library.

Depth selection

For depth, it's exactly the same as the time.

In [4]:
lon_min = -10
lon_max = 0
lat_min = 42
lat_max = 52
geotest = (datapgn.nav_lon > lon_min) & (datapgn.nav_lon < lon_max) & (datapgn.nav_lat > lat_min) & (datapgn.nav_lat < lat_max)
geoindex = np.argwhere(geotest.values)
xmin = min(geoindex[:,1])
xmax = max(geoindex[:,1])
ymin = min(geoindex[:,0])
ymax = max(geoindex[:,0])

In the cell below we extract data over our needed area, vertical levels and time.
We combine the use of the two methods presented before to select our data.
The slice command allows us to get all values between the 2 numbers.
For discrete extraction you can use a list (example below with deptht).
Note that xarray object use the dict formalism.

In [5]:
ldate = pd.date_range(start="20150301",end="20150307",freq="W-MON") # all monday between start and end
biscay = datapgn.isel({'x':slice(xmin,xmax),'y':slice(ymin,ymax)}).sel(deptht=[0,5,50,500],method='nearest').sel({'time_counter':ldate},method="nearest")

Dimensions of biscay dataset have been reduce to correct values.

In [6]:
biscay.dims
Out[6]:
Frozen(SortedKeysDict(OrderedDict([('y', 182), ('x', 127), ('deptht', 4), ('time_counter', 1)])))
In [7]:
biscay.nav_lon
Out[7]:
<xarray.DataArray 'nav_lon' (y: 182, x: 127)>
array([[-10.720358, -10.636135, -10.551914, ...,  -0.289115,  -0.205113,
         -0.121114],
       [-10.718001, -10.633769, -10.549537, ...,  -0.285598,  -0.201588,
         -0.117581],
       [-10.715626, -10.631382, -10.54714 , ...,  -0.282051,  -0.198033,
         -0.114018],
       ...,
       [ -9.904088,  -9.816092,  -9.728106, ...,   0.930388,   1.017018,
          1.103634],
       [ -9.896829,  -9.8088  ,  -9.720779, ...,   0.941264,   1.027917,
          1.114557],
       [ -9.889535,  -9.801472,  -9.713416, ...,   0.952194,   1.03887 ,
          1.125533]], dtype=float32)
Coordinates:
    nav_lon  (y, x) float32 -10.720358 -10.636135 ... 1.0388705 1.1255332
    nav_lat  (y, x) float32 42.050804 42.049053 42.047295 ... 51.846825 51.83977
  * x        (x) int32 3315 3316 3317 3318 3319 ... 3437 3438 3439 3440 3441
  * y        (y) int32 2052 2053 2054 2055 2056 ... 2229 2230 2231 2232 2233
Attributes:
    units:          degrees_east
    valid_min:      -179.99998963909385
    valid_max:      180.0
    long_name:      Longitude
    nav_model:      Default grid
    standard_name:  longitude

3 - Write file on disk (if needed)

In [8]:
biscay.to_netcdf("/homelocal/mydisk/Myfileonbiscay.nc")

4 - Crossing 73° longitude

The problem is that our previous test in order to get indexes from longitude and latitude will fails here.
It returns all indexes because 73° longitude is the first value of the array.
The change consists in the test on longitude.

In [9]:
lon_min = 65
lon_max = 75
lat_min = -54
lat_max = -44
geotest = ((datapgn.nav_lon < lon_min) | (datapgn.nav_lon > lon_max)) & (datapgn.nav_lat > lat_min) & (datapgn.nav_lat < lat_max)
geoindex = np.argwhere(geotest.values)
xmin = min(geoindex[:,1])
xmax = max(geoindex[:,1])
ymin = min(geoindex[:,0])
ymax = max(geoindex[:,0])
In [10]:
guelenker = datapgn.isel({'x':np.concatenate([np.arange(2,xmin+1),np.arange(xmax,datapgn.x.size)]),'y':slice(ymin,ymax)}).sel(deptht=[0,5,50,500],method='nearest').sel({'time_counter':ldate},method="nearest")

The dataset sort indexes by increasing order so we have the east of our area misplaced.
In order to correct this behaviour we use the roll method.

In [11]:
kerguelen = guelenker.roll(x=-xmin+1, roll_coords=True)

Second part : Get data with nco

We use nco after installing via conda install.

Note that you cannot extract by using longitude and latitude as they are variables and not dimensions.
You have to use index for dimension x and y so you have to know these indexes.
You can use both indexes or values for deptht and time_counter dimension as they are also variables.
For this test we select the 2 first depth levels by index -d deptht,0,1 and time_counter by values.
For help just write ncks.

Over Bay of Biscay

In [12]:
ncks -h -O -d x,3441,3442 -d y,2233,2234 -d deptht,0,1 -d time_counter,"2015-03-01 12:00:00","2015-03-02 12:00:00" -v nav_lon,nav_lat,deptht,votemper https://tds.mercator-ocean.fr/thredds/dodsC/psy4v3r1-daily-gridT biscay_nco.nc

Over Kerguelen (cross 73° longitude), you have to extract in two times with a first pass excluding the duplicate values and then extract across longitude indexes.

In [13]:
ncks -h -O -d x,2, -d y,721,725 -d deptht,0,1 -d time_counter,"2015-03-01 12:00:00","2015-03-02 12:00:00" -v nav_lon,nav_lat,deptht,votemper http://tds.mercator-ocean.fr/thredds/dodsC/psy4v3r1-daily-gridT guelenker_nco.nc
In [14]:
ncks -h -O -d x,4223,23 guelenker_nco.nc kerguelen_nco.nc