Contents

Using ADAM-API to access MODIS Aqua CHL

Contents

Using ADAM-API to access MODIS Aqua CHL#

you need to get an account to https://reliance.adamplatform.eu/ (use ORCID to authenticate) and key your ADAM API key
make sure you save your ADAM API key in a file $HOME/adam-key

!pip install adamapi

WARNING: The directory '/home/jovyan/.cache/pip' or its parent directory is not owned or is not writable by the current user. The cache has been disabled. Check the permissions and owner of that directory. If executing pip with sudo, you may want sudo's -H flag.
Requirement already satisfied: adamapi in /opt/conda/lib/python3.8/site-packages (2.0.11)
Requirement already satisfied: requests>=2.22.0 in /opt/conda/lib/python3.8/site-packages (from adamapi) (2.25.1)
Requirement already satisfied: imageio in /opt/conda/lib/python3.8/site-packages (from adamapi) (2.9.0)
Requirement already satisfied: certifi>=2017.4.17 in /opt/conda/lib/python3.8/site-packages (from requests>=2.22.0->adamapi) (2020.12.5)
Requirement already satisfied: chardet<5,>=3.0.2 in /opt/conda/lib/python3.8/site-packages (from requests>=2.22.0->adamapi) (4.0.0)
Requirement already satisfied: idna<3,>=2.5 in /opt/conda/lib/python3.8/site-packages (from requests>=2.22.0->adamapi) (2.10)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /opt/conda/lib/python3.8/site-packages (from requests>=2.22.0->adamapi) (1.26.4)
Requirement already satisfied: pillow in /opt/conda/lib/python3.8/site-packages (from imageio->adamapi) (8.1.2)
Requirement already satisfied: numpy in /opt/conda/lib/python3.8/site-packages (from imageio->adamapi) (1.20.1)

import os
import glob
import pathlib
import zipfile
import adamapi as adam
import xarray as xr
from datetime import datetime
import matplotlib
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cmaps

Authenticate to ADAM platform with ADAM API key#

adam_key = open(os.path.join(os.environ['HOME'],"adam-key")).read().rstrip()

a = adam.Auth()

a.setKey(adam_key)
a.setAdamCore('https://reliance.adamplatform.eu')
a.authorize()

{'expires_at': '2021-11-03T14:52:14.481Z',
 'access_token': '18d781c8642a45779a6a7d9df8b679fd',
 'refresh_token': '7ce4f620ce1848f9b34ae167d801de9b',
 'expires_in': 3600}

Discover MOD_Aqua datasets#

This step is useful to get the dataset identifier (unique for a given datacube)

def discoverDasasets(a, search_name):
    datasets = adam.Datasets(a)
    catalogue = datasets.getDatasets()
    #Extracting the size of the catalogue
    total = catalogue['properties']['totalResults']
    items = catalogue['properties']['itemsPerPage']
    pages = total//items
    
    print('----------------------------------------------------------------------')
    print('\033[1m' + 'List of available datasets:')
    print ('\033[0m')

    #Extracting the list of datasets across the whole catalogue
    for i in range(0,pages):
        page = datasets.getDatasets(page = i)
        for element in page['content']: 
            if search_name in element['title'] :
                print(element['title'] + "\033[1m" + " --> datasetId "+ "\033[0m" + "= " + element['datasetId'])
    return datasets

datasets = discoverDasasets(a, 'MOD_Aqua')

----------------------------------------------------------------------
List of available datasets:

Get metadata from Modis Aqua Chlorophylle Concentration#

datasetID = '69618:MODh20chlMO_4km'

print('\033[1;34m' + 'Metadata of ' + datasetID + ':')
print ('\033[0;0m')

paged = datasets.getDatasets(datasetID)
for i in paged.items():
    print("\033[1m" +  str(i[0]) + "\033[0m" + ': ' + str(i[1]))

Discover and select products from a dataset#

for a given time range and spatial coverage

Get data over Artic region#

The geometry field is extracted from a GeoJSON object , retrieving the value of the “feature” element.

Search data#

only print the first 10 products

from adamapi import Search

The GeoJson object needs to be rearranged according to the counterclockwise winding order.This operation is executed in the next few lines to obtain a geometry that meets the requirements of the method. Geom_1 is the final result to be used in the discovery operation.

geom_1 = "{'type': 'Polygon', 'coordinates': [[[-180, 90], [180, 90], [180, 60], [-180, 60], [-180, 90]]]}"

start_date = '2003-01-01'
end_date = '2003-12-31'

search = Search( a )
results = search.getProducts(
    datasetID, 
    geometry= geom_1,
    startDate = start_date,
    endDate = end_date
 )

# Printing the results

print('\033[1m' + 'List of available products:')
print ('\033[0m')
count = 1
for i in results['content']:

        print("\033[1;31;1m" + "#" + str(count))
        print ('\033[0m')
        for k in i.items():
            print(str(k[0]) + ': ' + str(k[1]))
        count = count+1
        print('------------------------------------')

Get data#

be aware that you alwasy get daily average from ADAM-API

def getZipData(auth, dataset_info):
    if not (pathlib.Path(pathlib.Path(dataset_info['outputFname']).stem).exists() or pathlib.Path(dataset_info['outputFname']).exists()):
        data = adam.GetData(auth)
        image = data.getData(
        datasetId = dataset_info['datasetID'],
        startDate = dataset_info['startDate'],
        endDate = dataset_info['endDate'],
        geometry = dataset_info['geometry'],
        outputFname = dataset_info['outputFname'])
        print(image)

%%time

output_file = './MOD_Aqua_mass_concentration_chlorophyll_concentration_in_sea_water_' + start_date + '-' + end_date + '.zip'

datasetInfo = {
    'datasetID'   : datasetID,
    'startDate'   : start_date,
    'endDate'     : end_date,
    'geometry'    : geom_1,
    'outputFname' : output_file
    }
getZipData(a, datasetInfo)

Data analysis and Visualization#

Unzip data#

def unzipData(filename):
    with zipfile.ZipFile(filename, 'r') as zip_ref:
        zip_ref.extractall(path = pathlib.Path(filename).stem)

if not pathlib.Path(pathlib.Path(output_file).stem).exists():
    unzipData(output_file)

Read data in xarray#

def paths_to_datetimeindex(paths):
    return  [datetime.strptime(date.split('_')[-1].split('.')[0], '%Y-%m-%dt%f') for date in paths]

def getData(dirtif, varname):
    geotiff_list = glob.glob(dirtif)
    # Create variable used for time axis
    time_var = xr.Variable('time', paths_to_datetimeindex(geotiff_list))
    # Load in and concatenate all individual GeoTIFFs
    geotiffs_da = xr.concat([xr.open_rasterio(i, parse_coordinates=True) for i in geotiff_list],
                        dim=time_var)
    # Covert our xarray.DataArray into a xarray.Dataset
    geotiffs_da = geotiffs_da.to_dataset('band')
    # Rename the dimensions to make it CF-convention compliant
    geotiffs_da = geotiffs_da.rename_dims({'y': 'latitude', 'x':'longitude'})
    # Rename the variable to a more useful name
    geotiffs_da = geotiffs_da.rename_vars({1: varname, 'y':'latitude', 'x':'longitude'})
    # set attribute to variable
    geotiffs_da[varname].attrs = {'units' : geotiffs_da.attrs[varname + '#units'], 'long_name' : geotiffs_da.attrs[varname + '#long_name']}
    return geotiffs_da

path_files = os.path.join(pathlib.Path(output_file).stem, '*.tif')
geotiff_ds = getData(path_files, 'chlor_a')

geotiff_ds

Visualization#

fig=plt.figure(figsize=(17,10))
# Define the projection
crs=ccrs.PlateCarree()

# We're using cartopy and are plotting in Orthographic projection 
# (see documentation on cartopy)
ax = plt.subplot(1, 1, 1, projection=ccrs.Mercator(central_longitude=12.0))
ax.coastlines(resolution='10m')

# We need to project our data to the new Mercator projection and for this we use `transform`.
# we set the original data projection in transform (here PlateCarree)
# we only plot values greather than 0
img = geotiff_ds['chlor_a'].isel(time=0).plot(ax=ax, transform=ccrs.PlateCarree(),  vmin=0, vmax=1, cmap=cmaps.BlueYellowRed)  

# Title for plot
plt.title('Mass concentration chlorophyll concentration in sea water \n',fontsize = 16, fontweight = 'bold', pad=10)
plt.savefig('Mod_Aqua-chlor_ARCTIC.png')

geotiff_ds['chlor_a'].to_netcdf('MOD_Aqua_Chl_arctic_2003.nc')