How to use the Mercator OpenDap
The access to this opendap requires authentication. To register please contact the Mercator-ocean service desk by :
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/.netrcAt the end, you must change the rights on the .netrc file using the following command:
chmod 600 .netrcFirst 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.
import numpy as np
import pandas as pd
import xarray as xr
datapgn = xr.open_dataset('https://tds.mercator-ocean.fr/thredds/dodsC/psy4v3r1-daily-gridT')
print(datapgn)
2 - Get dataset information : dimensions and variables¶
datapgn
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.
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.
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.
biscay.dims
biscay.nav_lon
3 - Write file on disk (if needed)¶
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.
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])
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.
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
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.
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
ncks -h -O -d x,4223,23 guelenker_nco.nc kerguelen_nco.nc