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')