Repository: OCHA-DAP/ocha-anticipy (language: python)

Issue:
Raster in top-level __init__
We should make the raster mixin importable from the top level and add to the readme, since we use it in a lot of pilot work
Specifying correct Python version when running pytest with pyenv
When testing the `pytest` framework, I ran into an issue where `pytest` didn't recognise the `--cov` argument passed to it. Below is a small example:

```shell
>> pytest --cov=src
ERROR: usage: pytest [options] [file_or_dir] [file_or_dir] [...]
pytest: error: unrecognized arguments: --cov=src
  inifile: None
  rootdir: /Users/caldwellst/Desktop/temp
```

Turns out users of `pyenv` could run into this issue where `pytest` runs off a different version of Python than what is active through `pyenv` (and that version in this case didn't have `pytest-cov` installed). Seems that there has been a fall through in my case with the shims, not sure why, so have to explicitly call `pytest` from my active Python version.

```shell
>> python -m pytest --cov=src
```
This is all based off of [this SO post](https://stackoverflow.com/questions/40718770/pytest-running-with-another-version-of-python). Raising this here in the case we wanted to somehow add this to documentation. It might be too niche to include, but wanted to at least have an issue active for others that may run into it.

Add message when data is overwritten (due to clobber)
Add to `check_file_existence` function something like:
`clobber_dict = {True: "", False: "not "}`
`f" Data already present in the folder will {clobber_dict[clobber]}be overwritten."`


Discussion / hints:

Fix patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 09341cec..81e73e9d 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,3 +1,6 @@
+default_language_version:
+    python: python3.11
+
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
     rev: v4.0.1
@@ -8,7 +11,7 @@ repos:
       - id: check-toml
       - id: check-yaml
   - repo: https://github.com/psf/black
-    rev: 21.9b0
+    rev: 22.3.0
     hooks:
       - id: black
   - repo: https://github.com/PyCQA/isort
@@ -36,31 +39,31 @@ repos:
           - types-PyYAML==5.4.10
           - types-setuptools==57.4.0
           - types-requests==2.25.9
+          - types-python-dateutil==2.8.19
         exclude: tests
   - repo: https://github.com/igorshubovych/markdownlint-cli
     rev: v0.29.0
     hooks:
       - id: markdownlint
   - repo: https://github.com/jazzband/pip-tools
-    rev: 6.3.0
+    rev: 6.9.0
     hooks:
       - id: pip-compile
         name: pip-compile requirements.txt
         files: setup.cfg
-        args: [setup.cfg, -q, --no-header,
-               -o, requirements/requirements.txt]
+        args: [setup.cfg, -q, -o, requirements/requirements.txt]
       - id: pip-compile
         name: pip-compile requirements-docs.txt
         files: requirements/docs.in
-        args: [ requirements/docs.in, -q, --no-header,
-                -o, requirements/requirements-docs.txt ]
+        args: [ requirements/docs.in, -q, -o,
+                requirements/requirements-docs.txt ]
       - id: pip-compile
         name: pip-compile requirements-tests.txt
         files: requirements/tests.in
-        args: [requirements/tests.in, -q, --no-header,
-               -o, requirements/requirements-tests.txt]
+        args: [requirements/tests.in, -q, -o,
+               requirements/requirements-tests.txt]
       - id: pip-compile
         name: pip-compile requirements-dev.txt
         files: requirements/(dev.in|docs.in|tests.in|requirements.txt)
-        args: [requirements/dev.in, -q, --no-header,
-               -o, requirements/requirements-dev.txt]
+        args: [requirements/dev.in, -q, -o,
+               requirements/requirements-dev.txt]
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 310dd3ef..7131157e 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -8,7 +8,7 @@ version: 2
 build:
   os: ubuntu-20.04
   tools:
-    python: "3.9"
+    python: "3.11"
 
 sphinx:
    configuration: docs/source/conf.py
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index f8e2d346..e5336f61 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -8,6 +8,36 @@ Changelog <https://keepachangelog.com/en/1.0.0/>`__, and this project
 adheres to `Semantic
 Versioning <https://semver.org/spec/v2.0.0.html>`__.
 
+[0.5.0] - 2022-12-02
+--------------------
+
+Added
+~~~~~
+
+- Modules for downloading and processing CHIRPS rainfall,
+  GloFAS river discharge, and USGS NDVI data
+- Utilities to streamline use of strings for dates across modules
+- COD AB configuration now has an ``admin{level}_name`` custom
+  layer name parameter
+- Logging for files being overwritten or not due to clobber
+
+Changed
+~~~~~~~
+
+- Documented and moved the raster processing module to the top level
+  for public access
+
+Removed
+~~~~~~~
+- Python 3.6 support
+
+Fixed
+~~~~~
+- The check in ``DataSource`` for the required configuration file
+  section now also checks if the section is ``None``
+- All available admin levels for DRC and Ethiopia are now accessible
+- IRI download method now checks request headers to verify authentication
+
 [0.4.2] - 2022-05-13
 --------------------
 
@@ -16,6 +46,7 @@ Fixed
 
 - Upgrade version of hdx-python-api to prevent bug when downloading
 
+
 [0.4.1] - 2022-05-10
 --------------------
 
diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
index 4b1927b5..afe62bf2 100644
--- a/CONTRIBUTING.rst
+++ b/CONTRIBUTING.rst
@@ -4,7 +4,15 @@ Development
 Environment
 -----------
 
-To set up the development environment, please install all packages from
+Development is currently done using Python 3.11. We recommend using a virtual
+environment such as ``venv``:
+
+.. code:: shell
+
+    python3.11 -m venv venv
+    source venv/bin/activate
+
+In your virtual environment, please install all packages from
 ``requirements/requirements-dev.txt``:
 
 .. code:: shell
@@ -33,14 +41,65 @@ To run the tests and view coverage, execute:
 
 .. code:: shell
 
-   pytest --cov=aatoolbox
+   python -m pytest --cov=aatoolbox
+
+Directly running pytest may
+[produce errors](https://stackoverflow.com/questions/40718770/pytest-running-with-another-version-of-python).
+by testing using a Python version or environment different from
+the local development environment, hence we recommend running
+as above.
 
 Documentation
 -------------
 
+Docstrings
+^^^^^^^^^^
+
 All public modules, classes and methods should be documented with
 `numpy-style <https://numpydoc.readthedocs.io/en/latest/format.html>`__
-docstrings, which can be rendered into HTML using the following command:
+docstrings.
+
+API documentation structure
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The API documentation structure is defined in ``api.rst``. The structure
+pulls from the docstrings for the modules, classes, and methods using
+Sphinx autodocumentation. ReST documentation is used to create navigable
+structure.
+
+Sphinx autodocumentation
+""""""""""""""""""""""""
+
+For a module and its children, the structure should follow the below guidelines.
+
+1. Document each module in the ``.py`` file itself. This is then documented
+directly using ``aautomodule``.
+
+2. Below this, document classes or functions as necessary. For classes
+that inherit directly from ``DataSource`` or other low-level base classes,
+just document ``:members:`` of the class. This ensures that only methods
+defined specifically for this class are documented.
+
+3. For classes that inherit from a base class that defines specific methods,
+like ``ChirpsDaily`` from ``_Chirps``, specify ``:inherited-members:`` as
+well. This ensures that only methods defined in the base class are documented.
+
+4. For something like the raster module classes, which inherits methods from
+a base class, but also defines unique methods within the subclass, use both
+``:members:`` and ``:inherited-members:``.
+
+Follow the example set out already in ``api.rst`` as you write the documentation.
+
+Module order
+""""""""""""
+
+The ordering of the ``api.rst`` documentation should match the ordering
+specified in the ``__all__`` list ``src/aatoolbox/__init__.py``.
+
+Build and view
+^^^^^^^^^^^^^^
+
+To build the documentation and test your implementation, use the following command:
 
 .. code:: shell
 
diff --git a/README.md b/README.md
index 555e6d1b..4f13e7d4 100644
--- a/README.md
+++ b/README.md
@@ -5,6 +5,7 @@
 [![PyPI Status](https://github.com/OCHA-DAP/pa-aa-toolbox/workflows/PyPI/badge.svg)](https://github.com/OCHA-DAP/pa-aa-toolbox/actions?query=workflow%3APyPI)
 [![Documentation Status](https://readthedocs.org/projects/aa-toolbox/badge/?version=latest)](https://aa-toolbox.readthedocs.io/en/latest/?badge=latest)
 [![Coverage Status](https://codecov.io/gh/OCHA-DAP/pa-aa-toolbox/branch/main/graph/badge.svg?token=JpWZc5js4y)](https://codecov.io/gh/OCHA-DAP/pa-aa-toolbox)
+[![pre-commit](https://img.shields.io/badge/pre--commit-enabled-brightgreen?logo=pre-commit)](https://github.com/pre-commit/pre-commit)
 [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black)
 [![Imports: isort](https://img.shields.io/badge/%20imports-isort-%231674b1?style=flat&labelColor=ef8336)](https://pycqa.github.io/isort/)
 
@@ -12,6 +13,15 @@ The [anticipatory action](https://centre.humdata.org/anticipatory-action/)
 (AA) toolbox is a Python package to support development of AA frameworks, by
 simplifying the downloading and processing of commonly used datasets.
 
+The datasets that we currently support are:
+
+- CHIRPS rainfall
+- COD ABs (Common Operational Datasets administrative boundaries)
+- FEWS NET food insecurity
+- GloFAS river discharge
+- IRI seasonal rainfall forecast
+- USGS NDVI (normalized difference vegetation index)
+
 For more information, please see the
 [documentation](https://aa-toolbox.readthedocs.io/en/latest/).
 
diff --git a/docs/source/api.rst b/docs/source/api.rst
index 644b469b..af5e2d1b 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -1,5 +1,142 @@
 AA Toolbox API
 ==============
 
-.. automodule:: aatoolbox
+Configuration
+-------------
+
+.. autoclass:: aatoolbox.CountryConfig
+
+.. autofunction:: aatoolbox.create_country_config
+.. autofunction:: aatoolbox.create_custom_country_config
+
+Data sources
+------------
+
+CHIRPS
+^^^^^^
+
+.. automodule:: aatoolbox.datasources.chirps.chirps
+
+Daily
+"""""
+
+.. autoclass:: aatoolbox.ChirpsDaily
+   :inherited-members:
+
+Monthly
+"""""""
+
+.. autoclass:: aatoolbox.ChirpsMonthly
+   :inherited-members:
+
+Common Operational Datasets
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: aatoolbox.datasources.codab.codab
+.. autoclass:: aatoolbox.CodAB
+   :members:
+
+FEWS NET
+^^^^^^^^
+
+.. automodule:: aatoolbox.datasources.fewsnet.fewsnet
+.. autoclass:: aatoolbox.FewsNet
+   :members:
+
+GloFAS
+^^^^^^
+
+.. automodule:: aatoolbox.datasources.glofas.glofas
+
+Reanalysis
+""""""""""
+
+.. autoclass:: aatoolbox.GlofasReanalysis
+   :inherited-members:
+
+Forecast
+""""""""
+
+.. autoclass:: aatoolbox.GlofasForecast
+   :inherited-members:
+
+Reforecast
+""""""""""
+
+.. autoclass:: aatoolbox.GlofasReforecast
+   :inherited-members:
+
+IRI
+^^^
+
+.. automodule:: aatoolbox.datasources.iri.iri_seasonal_forecast
+
+Probability forecast
+""""""""""""""""""""
+
+.. autoclass:: aatoolbox.IriForecastProb
+   :inherited-members:
+
+Dominant forecast
+"""""""""""""""""
+
+.. autoclass:: aatoolbox.IriForecastDominant
+   :inherited-members:
+
+NDVI (USGS eMODIS)
+^^^^^^^^^^^^^^^^^^
+
+.. automodule:: aatoolbox.datasources.usgs.ndvi_products
+
+Smoothed
+""""""""
+
+.. autoclass:: aatoolbox.UsgsNdviSmoothed
+   :inherited-members:
+
+Percent of median
+"""""""""""""""""
+
+.. autoclass:: aatoolbox.UsgsNdviPctMedian
+   :inherited-members:
+
+Median anomaly
+""""""""""""""
+
+.. autoclass:: aatoolbox.UsgsNdviMedianAnomaly
+   :inherited-members:
+
+Difference from previous year
+"""""""""""""""""""""""""""""
+
+.. autoclass:: aatoolbox.UsgsNdviYearDifference
+   :inherited-members:
+
+Utilities
+---------
+
+GeoboundingBox
+^^^^^^^^^^^^^^
+
+.. automodule:: aatoolbox.utils.geoboundingbox
+.. autoclass:: aatoolbox.GeoBoundingBox
+   :members:
+
+Raster module
+^^^^^^^^^^^^^
+
+.. automodule:: aatoolbox.utils.raster
+
+Data arrays
+"""""""""""
+
+.. autoclass:: aatoolbox.utils.raster.AatRasterArray
+   :members:
+   :inherited-members: RasterArray
+
+Datasets
+""""""""
+
+.. autoclass:: aatoolbox.utils.raster.AatRasterDataset
    :members:
+   :inherited-members: RasterDataset
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 2e7a128e..3775b14b 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -19,7 +19,7 @@
 
 project = "AA Toolbox"
 author = "Centre for Humanitarian Data Predictive Analytics Team"
-copyright = "2021, Centre for Humanitarian Data Predictive Analytics Team"
+copyright = "2022, Centre for Humanitarian Data Predictive Analytics Team"
 
 # -- General configuration ---------------------------------------------------
 
diff --git a/docs/source/config.rst b/docs/source/config.rst
index bf5aab24..84bd02f5 100644
--- a/docs/source/config.rst
+++ b/docs/source/config.rst
@@ -13,7 +13,9 @@ Built-in Country Configuration
 Several countries are already supported out-of-the box,
 and we are continuously adding more.
 
+.. _list of supported countries:
 A list of countries that we currently support:
+
 - Burkina Faso (BFA)
 - Bangladesh (BGD)
 - Democratic Republic of the Congo (COD)
diff --git a/docs/source/datasources.rst b/docs/source/datasources.rst
index e12523bc..1629a34c 100644
--- a/docs/source/datasources.rst
+++ b/docs/source/datasources.rst
@@ -6,9 +6,12 @@ AA Toolbox supports the following data sources:
 .. toctree::
    :maxdepth: 1
 
+   datasources/chirps
    datasources/codab
    datasources/fewsnet
+   datasources/glofas
    datasources/iri_seasonal_forecast
+   datasources/usgs_ndvi
 
 In general, each data source is available as a class with three methods:
 
diff --git a/docs/source/datasources/chirps.rst b/docs/source/datasources/chirps.rst
new file mode 100644
index 00000000..f7561832
--- /dev/null
+++ b/docs/source/datasources/chirps.rst
@@ -0,0 +1,123 @@
+CHIRPS
+======
+
+Background
+----------
+
+The Climate Hazards Center at UC Santa Barbara (`CHC <https://www.chc.ucsb.edu/>`_),
+in collaboration with scientists at the USGS Earth Resources Observation and Science
+(`EROS <https://www.usgs.gov/centers/eros>`_) Center, created the Climate Hazards group
+Infrared Precipitation with Stations (`CHIRPS <https://www.chc.ucsb.edu/data/chirps>`_)
+environmental record.
+
+CHIRPS is a quasi-global (50°S-50°N), high resolution (0.05°), daily, pentadal,
+and monthly precipitation dataset, ranging from 1981 to near-present. It incorporates
+CHPclim (CHC in-house climatology), 0.05° resolution satellite imagery, and in-situ station
+data to create gridded rainfall time series for trend analysis and seasonal drought
+monitoring. More information are available `here <https://www.nature.com/articles/sdata201566>`_.
+
+Through the AA toolbox, it is possible to have access to daily and monthly data.
+
+The data is downloaded from `IRI's maproom
+<http://iridl.ldeo.columbia.edu/SOURCES/.UCSB/.CHIRPS/.v2p0>`_ as the platform allows
+for the selection of geographical areas, and the delay in the upload
+of the data between CHC servers and IRI servers is less than one day.
+The data for a certain month is usually available two-three weeks after the end
+of the month.
+
+
+Usage
+-----
+
+In order to access CHIRPS data, you first need to create a country
+configuration for the country you would like to get data for. In this example
+we will use Burkina Faso, which corresponds to the ISO3 code "bfa":
+
+.. code-block:: python
+
+    from aatoolbox import create_country_config
+    country_config = create_country_config(iso3="bfa")
+
+Another required input is the geographic area of interest. A simple
+way to identify the area around the chosen country is to use the COD
+administrative boundaries. As a first step, you need to download them:
+
+.. code-block:: python
+
+    from aatoolbox import CodAB
+    codab = CodAB(country_config=country_config)
+    codab.download()
+    admin0 = codab.load(admin_level=0)
+
+Next, an instance of the class GeoBoundingBox needs to be created:
+
+.. code-block:: python
+
+    from aatoolbox import GeoBoundingBox
+    geo_bounding_box = GeoBoundingBox.from_shape(admin0)
+
+You can then choose start and end date of the dataset to be considered in your
+analysis, by specifying both of them as `datetime.date` objects.
+Both parameters are optional: if not specified, the start and end
+date will be respectively set to 1981-1-1 and the most recent date for which
+the data is available on the server.
+
+You can make use of the Chirps classes: two classes are
+available, one to retrieve daily data and one to retrieve monthly data. Daily
+data are available with two spatial resolutions (0.05 and 0.25 degrees),
+whereas monthly data can be obtained only with a 0.05-degree resolution.
+
+In the following example, monthly data ranging from February 2001 to March
+2006 is downloaded.
+
+.. code-block:: python
+
+    from aatoolbox import ChirpsMonthly
+
+    start_date = datetime.date(year=2001, month=2, day=1)
+    end_date = datetime.date(year=2006, month=3, day=31)
+
+    chirps_monthly = ChirpsMonthly(
+        country_config=country_config,
+        geo_bounding_box=geo_bounding_box,
+        start_date=start_date,
+        end_date=end_date
+        )
+
+    chirps_monthly.download()
+
+Similarly, the code snippet below allows to download daily CHIRPS data with
+0.25 degrees resolution and ranging from October 23, 2007 to the most recent
+available data:
+
+.. code-block:: python
+
+    from aatoolbox import ChirpsDaily
+
+    start_date = datetime.date(year=2007, month=10, day=23)
+
+    chirps_daily = ChirpsDaily(
+        country_config=country_config,
+        geo_bounding_box=geo_bounding_box,
+        start_date=start_date
+        )
+
+    chirps_daily.download()
+
+After having downloaded the datasets, a processing step is needed before
+being able to use them.
+
+.. code-block:: python
+
+    chirps_monthly.process()
+    chirps_daily.process()
+
+Finally, the data can be loaded as an ``xarray`` dataset, which is the result
+of the merging of all processed datasets, with fixed time resolution and
+location. Below are two examples of the use of the ``load`` method,
+respectively for monthly and daily data.
+
+.. code-block:: python
+
+    chirps_monthly_data = chirps_monthly.load()
+    chirps_daily_data = chirps_daily.load()
diff --git a/docs/source/datasources/codab.rst b/docs/source/datasources/codab.rst
index 5bbd5ecb..00726c2e 100644
--- a/docs/source/datasources/codab.rst
+++ b/docs/source/datasources/codab.rst
@@ -78,6 +78,7 @@ should be setup as follows:
       hdx_dataset_name: npl_admbnda_nd_20201117_SHP.zip
       layer_base_name: npl_admbnda_adm{admin_level}_nd_20201117.shp
       admin_level_max: 2
+      admin2_name: npl_admbnda_adm2_nd.shp
       custom_layer_names:
         - npl_admbnda_districts_nd_20201117
 
@@ -96,8 +97,13 @@ In the case of Nepal, the layers have the names ``npl_admbnda_adm0_nd_20201117.s
 ``npl_admbnda_adm1_nd_20201117.shp``, and ``npl_admbnda_adm2_nd_20201117.shp``
 
 ``admin_level_max``: The maximum admin level available in the layers. In the case of Nepal,
-the layer level numbers range from 0 to 2, so the maximum should be 2.
+the layer level numbers range from 0 to 2, so the maximum should be 2. In general the
+maximum admin level should not exceed 4.
 
-``custom_layer_name``: A place to list any other layers that don't correspond to the
+``admin{level}_name``: An optional parameter for any admin level (``level`` can range from 0 to 4)
+whose layer names do not match the ``layer_base_name`` pattern. This example for Nepal
+is contrived, but this issue does exist for COD ABs from countries such as Ethiopia and DRC.
+
+``custom_layer_name``: An optional place to list any other layers that don't correspond to the
 admin level format specified above. In the case of Nepal, there is a layer for districts
 with the name ``npl_admbnda_districts_nd_20201117``.
diff --git a/docs/source/datasources/glofas.rst b/docs/source/datasources/glofas.rst
new file mode 100644
index 00000000..317830bd
--- /dev/null
+++ b/docs/source/datasources/glofas.rst
@@ -0,0 +1,354 @@
+GloFAS
+======
+
+Background
+----------
+
+The
+`Global Flood Awareness System (GloFAS)
+<https://www.globalfloods.eu/>`_
+is part of the
+`Copernicus Emergency Management Service (CEMS).
+<https://emergency.copernicus.eu/>`_.
+It provides a global model and forecast of river flow
+for large catchments and riverine flooding,
+intended to be complementary to national hydrological and
+meteorological services.
+
+The processed and curated data,
+plus several other products, are available to view through the GloFAS
+`map viewer
+<https://www.globalfloods.eu/glofas-forecasting/>`_.
+While the map viewer makes image graphs of modelled and forecasted river
+discharge available, the actual data are needed to perform more in-depth
+analyses, such as computing specific return periods, knowing when
+a certain river discharge threshold was exceeded in the past, or
+detailed forecast performance at a specific location,
+
+The raw data in raster format can be downloaded from  the
+`Climate Data Store
+<https://cds.climate.copernicus.eu/#!/home>`_.
+However, the API can be tricky to navigate, with query size limits
+and long download times which depend on the query structure.
+Furthermore, the data are in the less commonly used GRIB
+file format, and the river discharge data at specific
+reporting points (those found on the map viewer)
+still need to be extracted. This GloFAS Python module takes
+care of all of these steps.
+
+A note about model versions: On 26 May 2021 GloFAS released `version 3
+<https://www.copernicus.eu/en/news/news/observer-whats-new-latest-glofas-31-release>`_
+of their model. While version 2 is still available on CDS, it
+will soon be phased out. Thus this module is currently only able
+to download version 3.
+
+Before describing the module usage, we will outline the different
+GloFAS datasets available: reanalysis, forecast, and reforecast.
+
+Reanalysis
+~~~~~~~~~~
+
+GloFAS-ERA5 reanalysis
+`(Harrigan et al.)
+<https://essd.copernicus.org/articles/12/2043/2020/>`_
+is a global gridded dataset of river discharge with
+a daily timestep and resolution of 0.1°,
+available from 1 January 1979 to the present day (updated daily).
+It is based on ECMWF's latest global atmospheric reanalysis
+`(ERA5, Herbasch et al.)
+<https://rmets.onlinelibrary.wiley.com/doi/10.1002/qj.3803>`_
+combined with the
+`Lisflood
+<https://ec-jrc.github.io/lisflood/>`_
+hydrological and channel routing model,
+and calibrated against observed river discharge
+`(Alfieri et al.)
+<https://www.sciencedirect.com/science/article/pii/S2589915519300331>`_.
+The data is available to download
+`here on CDS
+<https://cds.climate.copernicus.eu/cdsapp#!/dataset/cems-glofas-historical?tab=overview>`_.
+
+Forecast
+~~~~~~~~
+
+The GloFAS forecast
+`(Harrigan et al.)
+<https://hess.copernicus.org/preprints/hess-2020-532/>`_
+uses the above model, applied to the
+`ECMWF integrated forecast system (IFS) ensemble forecast
+<https://www.ecmwf.int/en/publications/ifs-documentation>`_.
+The IFS has 51 ensemble members out to a lead time of 15 days.
+The ensemble is composed of a single control member, generated from
+the most accurate estimate of current conditions, and 50 other
+members that have their initial conditions perturbed.
+
+Twice per week (Mondays and Thursdays), the IFS is extended to run up to 46 days ahead
+at a coarser resolution (~36 km). GloFAS makes available days
+16 to 30 for this extended forecast.
+
+Version 3 of the forecast data is available to download from 26 May 2021 until
+present day (updated daily) `here
+<https://cds.climate.copernicus.eu/cdsapp#!/dataset/cems-glofas-forecast?tab=overview>`_.
+While CDS does have version 3 pre-release data from 2020-2021,
+we understand that there were some small issues that were fixed
+in the operational version, so at this point in time this module
+only retrieves data from the operational model.
+
+Reforecast
+~~~~~~~~~~
+
+The GloFAS reforecast
+`(Harrigan et al.)
+<https://hess.copernicus.org/preprints/hess-2020-532/>`_
+is similar to the forecast, but uses the ECMWF IFS
+`reforecast
+<https://www.ecmwf.int/en/forecasts/documentation-and-support/extended-range/re-forecast-medium-and-extended-forecast-range>`_.
+The reforecast is initialized twice per week (Monday and Thursdays)
+and has 11 ensemble members.
+The data runs from January 1999 to December 2018,
+and is available to download
+`from CDS here
+<https://cds.climate.copernicus.eu/cdsapp#!/dataset/cems-glofas-reforecast?tab=overview>`_.
+
+Usage
+-----
+
+This module performs three basic steps:
+
+#. Download the raw GloFAS data, in GRIB file format.
+   Each file is a raster based on a
+   user-specified region of interest, and contains data for a full
+   day, month, or year, depending on the GloFAS product.
+#. Process the data. For each raw GRIB file,
+   river discharge for a set of user-provided pixel locations is extracted
+   and stored in NetCDF format.
+#. Load the data. The processed files, which are split by time in the same
+   way as the raw data,
+   are read in as a single `xarray.DataSet`.
+
+CDS API
+~~~~~~~
+
+The download step makes use of the
+`CDS API
+<https://cds.climate.copernicus.eu/api-how-to>`_.
+You will need to register an account on CDS, then once logged in go to your
+`user page
+<https://cds.climate.copernicus.eu/user/>`_
+and note down your UID and API key.
+
+The
+`how-to page
+<https://cds.climate.copernicus.eu/api-how-to>`_
+has instructions for each operating system, but we summarize the instructions below:
+
+Windows: Create a file called `%USERPROFILE%\.cdsapirc` where, `%USERPROFILE%` is usually located
+in your `C:\Users\Username` folder.
+
+Mac and Linux: Create the file `~/.cdsapirc`.
+
+Then, for all operating systems, add the following to the created file:
+
+.. code-block:: shell
+
+    url: https://cds.climate.copernicus.eu/api/v2
+    key: ${UID}:${API-key}
+
+where `$UID` is your UID and `$API-key` is your API key.
+
+ecCodes
+~~~~~~~
+
+To read in .GRIB files using `xarray`, you'll need to install
+`ecCodes
+<https://confluence.ecmwf.int/display/ECC/What+is+ecCodes>`_
+on your machine.
+For all operating systems, this can be done through
+the
+`source distribution
+<https://confluence.ecmwf.int/display/ECC/ecCodes+installation>`_,
+or in a `conda` environment with the command:
+
+.. code-block:: shell
+
+    conda install -c conda-forge eccodes
+
+Some OS-specific binaries are also available. For Linux,
+`python3-eccodes` can be found as a `.deb` and `.rpm`
+(check your specific Linux distribution for the latest version).
+For Mac, according to the
+`ecCodes-Python documentation
+<https://github.com/ecmwf/eccodes-python#system-dependencies>_`,
+ecCodes can be installed using `brew`:
+
+.. code-block:: shell
+
+    brew install eccodes
+
+Reporting points
+~~~~~~~~~~~~~~~~
+
+Next, if it :ref:`doesn't already exist<list of supported countries>`,
+you need to create a country configuration
+for the country you would like to analyze.
+
+An example country config for Bangladesh is:
+
+.. code-block:: yaml
+
+    iso3: bgd
+    glofas:
+      reporting_points:
+      - name: Bahadurabad
+        lon: 89.65
+        lat: 25.15
+      - name: Hardinge Bridge
+        lon: 89.05
+        lat: 24.05
+
+The reporting points indicate the raster file coordinates used
+to extract the river discharge for a particular location.
+Those in the above example have been taken from the
+`GloFAS map viewer
+<https://www.globalfloods.eu/glofas-forecasting/>`_.
+If you select "Reporting Points" from the "Hydrological"
+menu at the top, they will appear as dots on the map. If you then
+click on one of the points, you are able to see
+information such as the station name, and LISFLOOD X and Y, which are
+the respective longitude and latitude used in the configuration file.
+
+Reporting point coordinates are manually selected by the GloFAS team to
+be representative of physical gauge locations, and to be located on a river
+in the model raster file. In principle, one could
+specify any set of coordinates that exists on the raster, but caution is advised
+when doing so.
+
+
+Running the code
+~~~~~~~~~~~~~~~~
+
+You can initialize a built-in country config as follows:
+
+.. code-block:: python
+
+    from aatoolbox import create_country_config
+
+    country_config = create_country_config(iso3="bgd")
+
+Another required input is the geographic area of interest, which will
+define the bounds of raw raster data to be downloaded. A simple
+way to identify the area around the chosen country is to use the COD
+administrative boundaries. You will need to download the data,
+and extract a geo bounding box:
+
+.. code-block:: python
+
+    from aatoolbox import CodAB, GeoBoundingBox
+
+    codab = CodAB(country_config=country_config)
+    codab.download()
+    admin0 = codab.load()
+    geo_bounding_box = GeoBoundingBox.from_shape(admin0)
+
+Note that the reporting points in the configuration file need to lie within
+the geographic area of interest.
+
+Next you need to instantiate the GloFAS class with the country config. For this
+example, we will use the GloFAS forecast, however the steps are nearly identical for the
+the reanalysis and reforecast -- the only differences are the acceptable date ranges,
+and that `leadtime_max` is **not** an input parameter to the reanlaysis.
+
+In this case, we would like the data for the past month. In general, we suggest
+specifying dates explicitly, as using e.g. `date.today()` will not produce
+the same results when run on a different day.
+
+.. code-block:: python
+
+    from datetime import date
+
+    from aatoolbox import GlofasForecast
+
+    glofas_forecast = GlofasForecast(
+        country_config=country_config,
+        geo_bounding_box=geo_bounding_box,
+        leadtime_max=15,
+        start_date=date(year=2022, month=9, day=22)
+        end_date=date(year=2022, month=10, day=22),
+     )
+
+We then need to download the GloFAS data. The module will download all the data
+between `start_date` and `end_date` (inclusive). The raw files are in .GRIB format,
+and are separated by day (forecast), month (reforecast), or year (reanalysis)
+depending on the data type, due to CDS query size limits.
+
+.. code-block:: python
+
+    glofas_forecast.download()
+
+The downloading process works by generating all requests required for the
+specified timeframe, sending them to CDS, and saving the request numbers in memory.
+You can see all your requests on the
+`CDS website
+<https://cds.climate.copernicus.eu/cdsapp#!/yourrequests>`_,
+and cancel any pending ones in case the run is interrupted (we hope to implement
+request number caching in a later version).
+
+The module then pings the CDS API every minute to check which requests have
+completed, and downloads those that have. This continues until all requests
+have been downloaded.
+
+This process can unfortunately take a long time, however, with this module
+we've tried to optimize the queries to be as fast as possible. In our experience,
+downloading the full reanalysis (from 1979 to today) is the fastest and takes a
+couple of hours, while the full reforecast (1999 to 2018) is the slowest and takes
+around a day to complete. It also depends how busy the queue is, which
+you can check
+`here
+<https://cds.climate.copernicus.eu/live/queue>`_.
+
+The next step is to process the files. First, we want to convert from the
+less used GRIB format to the more common and flexible NetCDF. Furthermore, rather
+than having the full raster, we extract the river discharge data at the
+reporting poitns specified in the configuration file. This can be done in a single
+step:
+
+.. code-block:: python
+
+    glofas_forecast.process()
+
+Note that each individual raw GRIB file is converted to a corresponding
+processed NetCDF file. This is to simplify the downloading and processing of
+addition data, i.e. for adding new dates.
+
+To load all of the processed files into a single dataframe (which can then
+be saved to a single NetCDF file or other compatible format), execute:
+
+.. code-block:: python
+
+    bgd_glofas_forecast_reporting_points = glofas_forecast.load()
+
+The full codde snippet is below:
+
+.. code-block:: python
+
+    from datetime import date
+
+    from aatoolbox import create_country_config, CodAB, \
+        GeoBoundingBox, GlofasForecast
+
+    codab = CodAB(country_config=country_config)
+    codab.download()
+    admin0 = codab.load()
+    geo_bounding_box = GeoBoundingBox.from_shape(admin0)
+
+    glofas_forecast = GlofasForecast(
+        country_config=country_config,
+        geo_bounding_box=geo_bounding_box,
+        leadtime_max=15,
+        start_date=date(year=2022, month=9, day=22)
+        end_date=date(year=2022, month=10, day=22),
+     )
+    glofas_forecast.download()
+    glofas_forecast.process()
+
+    bgd_glofas_forecast_reporting_points = glofas_forecast.load()
diff --git a/docs/source/datasources/usgs_ndvi.rst b/docs/source/datasources/usgs_ndvi.rst
new file mode 100644
index 00000000..1bcf98b0
--- /dev/null
+++ b/docs/source/datasources/usgs_ndvi.rst
@@ -0,0 +1,142 @@
+NDVI (USGS eMODIS)
+==================
+
+Warning
+-------
+
+The MODIS sensor has been reported by USGS to
+have degraded in quality since May 2022 (dekad 13), and
+updates to this data source have stopped. This module
+remains for users to have access to historic data but
+recent data is unavailable and care should be used
+analyzing any data since dekad 13 of 2022.
+
+Background
+----------
+
+The United States Geological Survey, `USGS <https://www.usgs.gov/>`_ provides a wide
+range of geo-spatial data and satellite imagery products in support of the
+Famine Early Warning Systems Network, `FEWS NET <https://earlywarning.usgs.gov/fews>`_.
+These products are used to support FEWS NET drought monitoring efforts, from rainfall
+data to vegetation measurements.
+
+The availability of data for each source varies,
+but are often available for every dekad or pentad for a number of years. The
+geographic coverage is typically comprehensive of the African continent, where
+the majority of FEWS NET's work is focused, with additional areas covered on an
+ad hoc basis. This module for now just provides access to the normalized
+difference vegetation index (`NDVI <https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index>`_)
+generated from eMODIS AQUA provided by the USGS.
+
+NDVI data is generated from eMODIS AQUA and published data includes temporally smoothed NDVI, median anomaly,
+difference from the previous year, and median anomaly presented as a percentile. Full methodological
+details are available on the `Documentation page <https://earlywarning.usgs.gov/fews/product/449#documentation>`_
+for the specific coverage area. The data is made available for the following areas of coverage:
+
+- `north-africa <https://earlywarning.usgs.gov/fews/product/449>`_
+- `east-africa <https://earlywarning.usgs.gov/fews/product/448>`_
+- `southern-africa <https://earlywarning.usgs.gov/fews/product/450>`_
+- `west-africa <https://earlywarning.usgs.gov/fews/product/451>`_
+- `central-asia <https://earlywarning.usgs.gov/fews/product/493>`_
+- `yemen <https://earlywarning.usgs.gov/fews/product/502>`_
+- `central-america <https://earlywarning.usgs.gov/fews/product/445>`_
+- `hispaniola <https://earlywarning.usgs.gov/fews/product/446>`_
+
+Data is published at the dekadal level, and is released soon after the end of the dekad.
+After a period of 3 dekads, data is updated with temporal smoothing and error correction
+for cloud cover. Files for a specific dekad and region can range from 30MB up to over 100MB,
+so downloading and processing can take a long time.
+
+
+USGS publishes its data through its `web portal <https://earlywarning.usgs.gov/fews/datadownloads/East%20Africa/eMODIS%20NDVI%20C6>`_,
+which allows downloading data for a single geographical area and dekad. This data
+is extracted from the `back end data explorer <https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/fews/web/africa/east/dekadal/emodis/ndvi_c6/temporallysmoothedndvi/downloads/monthly/>`_.
+This module is designed to allow programmatic access to and analysis of NDVI data from the data explorer.
+
+Usage
+-----
+
+To use this class, you first need to create a country configuration
+for the country you would like to analyze. You also need to supply the area
+of coverage that the country data is contained in. You can Make sure FEWS NET covers the country
+of interest. You can find area of coverage for NDVI by searching for eMODIS NDVI C6 on
+the  `USGS search portal <https://earlywarning.usgs.gov/fews/search>`_.
+The valid values of the area names are listed above under Background.
+
+An example country config for Ethiopia is:
+
+.. code-block:: python
+
+    iso3: eth
+    usgs_ndvi:
+        region_name: east-africa
+
+For a number of countries this config has been implemented already.
+For others, you can create a custom country config.
+Once you have created the config you can load it:
+
+.. code-block:: python
+
+    from aatoolbox import create_country_config
+    country_config = create_country_config(iso3="eth")
+
+Next you need to instantiate the USGS NDVI class with the country config. For this
+example, we will use the class for accessing temporally smoothed NDVI values. Since
+the size of the data we will be working with is so large, we will only look at a
+sample of data from a few dekads, corresponding to the dates in January:
+
+.. code-block:: python
+
+    from aatoolbox import UsgsNdviSmoothed
+    ndvi_smooth = UsgsNdviSmoothed(
+        country_config=country_config,
+        start_date="2021-01-01",
+        end_date="2021-01-31"
+    )
+
+We will first need to download the NDVI data if you haven't used the class before.
+The code downloads all data between `start_date` and `end_date`. Be careful if not
+setting start and end date during instantiation, the method will automatically
+look for data between the first available data in the 19th dekad of 2002 and today.
+We simply have to call the method to download available data.
+
+.. code-block:: python
+
+    ndvi_smooth.download()
+
+Since each NDVI file is downloaded as a .tif, we can load them in individually as
+rasters. We can do this for a specific date. This time we will specify the date
+as a year and dekad tuple, although specify as a `datetime.date` or string as
+above is equally appropriate.
+
+.. code-block:: python
+
+    ndvi_2021_01 = ndvi_smooth.load_raster(date=(2021, 1))
+
+If time series analysis is desired on the NDVI data, the user can manually do this
+using loaded in data arrays. However, given the size of the individual raster files,
+this module provides the `process()` method that calculate statistics for
+a given area. For this, we need to provide a geodataframe and column to aggregate to.
+
+Let's load the country administrative boundaries for Ethiopia from our country
+config, and use this to calculate basic statistics for the dekads we've loaded.
+
+.. code-block:: python
+
+    from aatoolbox import CodAB
+
+    eth_cod = CodAB(country_config=country_config)
+    # assuming you've downloaded the file already
+    eth_gdf0 = eth_cod.load(admin_level=0)
+    ndvi_smooth.process(
+        gdf=eth_gdf0,
+        feature_col="ADM0_EN"
+    )
+
+We can then load in and use the calculated statistics using `load()`. Since the
+processed file is saved based on the `feature_col` name, we have to pass
+`feature_col` to load in the correct data.
+
+.. code-block:: python
+
+    ndvi_smooth.load(feature_col="ADM0_EN")
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 62378cc9..62de243a 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -28,6 +28,7 @@ additional functionality.
    quickstart
    config
    datasources
+   utilities
 
 API Reference
 -------------
diff --git a/docs/source/introduction.rst b/docs/source/introduction.rst
index 12ac5518..39d0feb3 100644
--- a/docs/source/introduction.rst
+++ b/docs/source/introduction.rst
@@ -12,6 +12,15 @@ This is where the AA Toolbox comes in. The goal of the toolbox is to to
 support development of AA frameworks by simplifying the downloading and
 processing of commonly used datasets.
 
+Datasets currently supported are:
+
+- CHIRPS rainfall
+- COD ABs (Common Operational Datasets administrative boundaries)
+- FEWS NET food insecurity
+- GloFAS river discharge
+- IRI seasonal rainfall forecast
+- USGS NDVI (normalized difference vegetation index)
+
 At the `Centre for Humanitarian Data <https://centre.humdata.org/>`_,
 we've now used AA Toolbox to develop
 some of our own AA frameworks. See the links below for more info:
diff --git a/docs/source/quickstart.rst b/docs/source/quickstart.rst
index 56a83658..a5c8e923 100644
--- a/docs/source/quickstart.rst
+++ b/docs/source/quickstart.rst
@@ -17,7 +17,7 @@ like the data to go.
 Install AA Toolbox
 ^^^^^^^^^^^^^^^^^^
 
-AA Toolbox supports Python 3.6 and newer. It can be installed in your Python
+AA Toolbox supports Python 3.7 and newer. It can be installed in your Python
 environment using the following command:
 
 .. code-block:: sh
diff --git a/docs/source/utilities.rst b/docs/source/utilities.rst
new file mode 100644
index 00000000..ad2132a9
--- /dev/null
+++ b/docs/source/utilities.rst
@@ -0,0 +1,9 @@
+Utilities
+=========
+
+AA Toolbox provides the following utilities:
+
+.. toctree::
+   :maxdepth: 1
+
+   utilities/raster
diff --git a/docs/source/utilities/raster.rst b/docs/source/utilities/raster.rst
new file mode 100644
index 00000000..72271a87
--- /dev/null
+++ b/docs/source/utilities/raster.rst
@@ -0,0 +1,160 @@
+Raster
+======
+
+Background
+----------
+
+This module extends xarray and rioxarray to provide
+additional functionality for raster processing and
+post-processing. These are used throughout the
+datasources to process various raster data, but are
+also useful directly for the end user, and is thus
+available for direct use. The extension is based on the
+guidance for how to `extend xarray
+<http://xarray.pydata.org/en/stable/internals/extending-xarray.html>`_.
+
+Usage
+-----
+
+Accessor
+^^^^^^^^
+
+To make the extension accessible, you only need to import
+``aatoolbox`` directly. Then, extension methods and properties
+are simply accessed using ``da.aat.method()`` or
+``da.aat.property`` on a data array or dataset.  As a simple
+example, we can create a data array that has
+incorrect coordinates (ascending latitude and descending
+longitude). Then we can use the extension to invert them to
+the correct ordering.
+
+.. code-block:: python
+
+    import xarray
+    import numpy
+    import aatoolbox
+
+    da = xarray.DataArray(
+      numpy.arange(16).reshape(4,4),
+      coords={"lat":numpy.array([87, 88, 89, 90]),
+              "lon":numpy.array([70, 69, 68, 67])}
+    )
+    da_inv = da.aat.invert_coordinates()
+
+    # check they have inverted
+    da_inv.get_index("lon")
+    #> Int64Index([67, 68, 69, 70], dtype='int64', name='lon')
+    da_inv.get_index("lat")
+    #> Int64Index([90, 89, 88, 87], dtype='int64', name='lat')
+
+Raster statistics
+^^^^^^^^^^^^^^^^^
+
+One of the more useful functionalities of the module
+is providing an easy ability to calcuate raster statistics
+across a specified set of geometries defined in a geodataframe.
+Paired with the administrative boundaries and raster data sources
+available through this library, we can easily calculate
+raster statistics.
+
+A full snippet of example code is available below.
+
+.. code-block:: python
+
+    import aatoolbox
+    import datetime
+
+    # load the administrative boundaries
+    country_config = aatoolbox.create_country_config(iso3="eth")
+    codab = aatoolbox.CodAB(country_config)
+    codab.download()
+    codab_eth = codab.load(admin_level=2)
+
+    # get geobounding box for CHIRPS downloading
+    geo_bounding_box = aatoolbox.GeoBoundingBox.from_shape(codab_eth)
+
+    # load CHIRPs data for processing
+    start_date = datetime.date(year=2001, month=2, day=1)
+    end_date = datetime.date(year=2006, month=3, day=31)
+
+    chirps_monthly = ChirpsMonthly(
+        country_config=country_config,
+        geo_bounding_box=geo_bounding_box,
+        start_date=start_date,
+        end_date=end_date
+        )
+
+    chirps_monthly.download()
+    chirps_monthly.process()
+    chirps_monthly_data = chirps_monthly.load()
+
+    # compute raster statistics
+
+    chirps_monthly_data.aat.compute_raster_stats(
+      gdf=codab_eth,
+      feature_col="ADM2_PCODE"
+    )
+
+Properties
+----------
+
+The user should be careful when accessing attributes of
+a data array when using the raster module. This module
+builds on `rioxarray <https://corteva.github.io/rioxarray>_`
+extensions, and thus methods and attributes accessible
+via ``da.rio.method()``  or ``da.rio.property`` are
+also accessible using ``da.aat.method()`` or
+``da.aat.property``. However, original rioxarray properties
+should be accessed using ``da.rio.property``.
+
+Let's create a simple data array where we want to specify
+the spatial dimensions explicitly because the coordinate
+names are not automatically detected.
+
+.. code-block:: python
+
+    import xarray
+    import numpy
+    import aatoolbox
+
+    da = xarray.DataArray(
+        numpy.arange(16).reshape(4,4),
+        coords={"a":numpy.array([90, 89, 88, 87]),
+                "b":numpy.array([70, 69, 68, 67])}
+    )
+
+We can set the spatial dimensions using
+``da.rio.set_spatial_dims()`` or call it directly
+from ``da.aat``.
+
+.. code-block:: python
+
+  da_new = da.aat.set_spatial_dims(
+    x_dim="a",
+    y_dim="b"
+  )
+
+However, even though we can set the dimensions
+using either accessor, we have to be careful
+accessing the properties.
+
+.. code-block:: python
+
+  da_new.rio.x_dim
+  #> 'a'
+
+  da_new.aat.x_dim
+  #> MissingSpatialDimensionError: x dimension not found.
+  #> 'rio.set_spatial_dims()' or using 'rename()' to change
+  #> the dimension name to 'x' can address this.
+
+Even though the method was called using ``aat``, the property
+is not accessible through it. Users need to be careful about
+accessing rioxarray properties using the ``aat`` accessor.
+
+For best practice, rioxarray methods and properties should all
+be accessed using ``rio``. These properties are ``rio.x_dim``,
+``rio.y_dim``, ``rio.shape``, ``rio.width``, ``rio.height``, and
+``rio.crs``. This module's methods and properties should be
+accessed using the ``aat`` accessor. These properties are
+``aat.t_dim`` and ``aat.longitude_range``.
diff --git a/requirements/requirements-dev.txt b/requirements/requirements-dev.txt
index dfd891d1..470251e3 100644
--- a/requirements/requirements-dev.txt
+++ b/requirements/requirements-dev.txt
@@ -1,3 +1,9 @@
+#
+# This file is autogenerated by pip-compile with python 3.11
+# To update, run:
+#
+#    pip-compile --output-file=requirements/requirements-dev.txt requirements/dev.in
+#
 affine==2.3.1
     # via
     #   -r requirements/requirements.txt
@@ -10,34 +16,36 @@ appdirs==1.4.4
     # via
     #   -r requirements/requirements.txt
     #   requests-cache
-attrs==21.4.0
+attrs==22.1.0
     # via
     #   -r requirements/requirements-tests.txt
     #   -r requirements/requirements.txt
     #   cattrs
+    #   cfgrib
+    #   eccodes
     #   fiona
     #   jsonlines
     #   jsonschema
     #   pytest
     #   rasterio
     #   requests-cache
-babel==2.10.1
+babel==2.11.0
     # via
     #   -r requirements/requirements-docs.txt
     #   sphinx
-basicauth==0.4.1
-    # via
-    #   -r requirements/requirements.txt
-    #   hdx-python-utilities
-bokeh==2.4.2
+bokeh==2.4.3
     # via
     #   -r requirements/requirements.txt
     #   dask
-cattrs==1.10.0
+build==0.9.0
+    # via pip-tools
+cattrs==22.2.0
     # via
     #   -r requirements/requirements.txt
     #   requests-cache
-certifi==2021.10.8
+cdsapi==0.5.1
+    # via -r requirements/requirements.txt
+certifi==2022.9.24
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements.txt
@@ -45,21 +53,25 @@ certifi==2021.10.8
     #   pyproj
     #   rasterio
     #   requests
-cffi==1.15.0
+cffi==1.15.1
     # via
     #   -r requirements/requirements.txt
     #   cryptography
+    #   eccodes
+cfgrib==0.9.10.3
+    # via -r requirements/requirements.txt
 cfgv==3.3.1
     # via pre-commit
-cftime==1.6.0
+cftime==1.6.2
     # via
     #   -r requirements/requirements-tests.txt
+    #   -r requirements/requirements.txt
     #   netcdf4
-chardet==4.0.0
+chardet==5.0.0
     # via
     #   -r requirements/requirements.txt
     #   frictionless
-charset-normalizer==2.0.12
+charset-normalizer==2.1.1
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements.txt
@@ -68,11 +80,13 @@ ckanapi==4.7
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-api
-click==8.1.2
+click==8.1.3
     # via
     #   -r requirements/requirements.txt
+    #   cfgrib
     #   click-plugins
     #   cligj
+    #   dask
     #   distributed
     #   fiona
     #   pip-tools
@@ -88,24 +102,28 @@ cligj==0.7.2
     #   -r requirements/requirements.txt
     #   fiona
     #   rasterio
-cloudpickle==2.0.0
+cloudpickle==2.2.0
     # via
     #   -r requirements/requirements.txt
     #   dask
     #   distributed
-colorama==0.4.4
+colorama==0.4.6
     # via
     #   -r requirements/requirements.txt
     #   typer
-coverage[toml]==6.3.2
+commonmark==0.9.1
+    # via
+    #   -r requirements/requirements.txt
+    #   rich
+coverage[toml]==6.5.0
     # via
     #   -r requirements/requirements-tests.txt
     #   pytest-cov
-cryptography==36.0.2
+cryptography==38.0.3
     # via
     #   -r requirements/requirements.txt
     #   pyopenssl
-dask==2022.04.1
+dask[complete]==2022.11.1
     # via
     #   -r requirements/requirements.txt
     #   distributed
@@ -115,15 +133,15 @@ decorator==5.1.1
     #   -r requirements/requirements.txt
     #   jsonpath-ng
     #   validators
-defopt==6.3.0
+defopt==6.4.0
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-api
-distlib==0.3.4
+distlib==0.3.6
     # via
     #   -r requirements/requirements-tests.txt
     #   virtualenv
-distributed==2022.04.1
+distributed==2022.11.1
     # via
     #   -r requirements/requirements.txt
     #   dask
@@ -143,7 +161,11 @@ docutils==0.17.1
     #   defopt
     #   sphinx
     #   sphinx-rtd-theme
-email-validator==1.1.3
+eccodes==1.5.0
+    # via
+    #   -r requirements/requirements.txt
+    #   cfgrib
+email-validator==1.3.0
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-api
@@ -151,36 +173,36 @@ et-xmlfile==1.1.0
     # via
     #   -r requirements/requirements.txt
     #   openpyxl
-exchangerates==0.3.4
-    # via
-    #   -r requirements/requirements.txt
-    #   hdx-python-country
-filelock==3.6.0
+filelock==3.8.0
     # via
     #   -r requirements/requirements-tests.txt
     #   tox
     #   virtualenv
-fiona==1.8.21
+findlibs==0.0.2
+    # via
+    #   -r requirements/requirements.txt
+    #   eccodes
+fiona==1.8.22
     # via
     #   -r requirements/requirements.txt
     #   geopandas
-frictionless[excel,json]==4.34.0
+frictionless==4.40.8
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-utilities
-fsspec==2022.3.0
+fsspec==2022.11.0
     # via
     #   -r requirements/requirements.txt
     #   dask
-geopandas==0.10.2
+geopandas==0.12.1
     # via -r requirements/requirements.txt
-hdx-python-api==5.6.5
+hdx-python-api==5.9.2
     # via -r requirements/requirements.txt
-hdx-python-country==3.2.1
+hdx-python-country==3.4.3
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-api
-hdx-python-utilities==3.3.1
+hdx-python-utilities==3.4.8
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-country
@@ -188,9 +210,9 @@ heapdict==1.0.1
     # via
     #   -r requirements/requirements.txt
     #   zict
-identify==2.4.12
+identify==2.5.9
     # via pre-commit
-idna==3.3
+idna==3.4
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements.txt
@@ -199,20 +221,12 @@ idna==3.3
 ijson==3.1.4
     # via
     #   -r requirements/requirements.txt
-    #   frictionless
-imagesize==1.3.0
-    # via
-    #   -r requirements/requirements-docs.txt
-    #   sphinx
-importlib-metadata==4.11.3
+    #   hdx-python-utilities
+imagesize==1.4.1
     # via
     #   -r requirements/requirements-docs.txt
     #   sphinx
-importlib-resources==5.7.1
-    # via
-    #   -r requirements/requirements.txt
-    #   jsonschema
-inflect==5.5.2
+inflect==6.0.2
     # via
     #   -r requirements/requirements.txt
     #   quantulum3
@@ -224,7 +238,7 @@ isodate==0.6.1
     # via
     #   -r requirements/requirements.txt
     #   frictionless
-jinja2==3.1.1
+jinja2==3.1.2
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements.txt
@@ -233,36 +247,37 @@ jinja2==3.1.1
     #   distributed
     #   frictionless
     #   sphinx
-jsonlines==3.0.0
+jsonlines==3.1.0
     # via
     #   -r requirements/requirements.txt
-    #   frictionless
+    #   hdx-python-utilities
 jsonpath-ng==1.5.3
     # via
     #   -r requirements/requirements.txt
     #   libhxl
-jsonschema==4.4.0
+jsonschema==4.17.1
     # via
     #   -r requirements/requirements.txt
     #   frictionless
     #   tableschema-to-template
-libhxl==4.25
+libhxl==4.27
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-country
 locket==1.0.0
     # via
     #   -r requirements/requirements.txt
+    #   distributed
     #   partd
 loguru==0.6.0
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-utilities
-lxml==4.8.0
+makefun==1.15.0
     # via
     #   -r requirements/requirements.txt
-    #   exchangerates
-marko==1.2.0
+    #   hdx-python-api
+marko==1.2.2
     # via
     #   -r requirements/requirements.txt
     #   frictionless
@@ -271,7 +286,7 @@ markupsafe==2.1.1
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements.txt
     #   jinja2
-msgpack==1.0.3
+msgpack==1.0.4
     # via
     #   -r requirements/requirements.txt
     #   distributed
@@ -279,78 +294,80 @@ munch==2.5.0
     # via
     #   -r requirements/requirements.txt
     #   fiona
-mypy-extensions==0.4.3
-    # via
-    #   -r requirements/requirements.txt
-    #   typing-inspect
 ndg-httpsclient==0.5.1
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-api
-netcdf4==1.5.8
-    # via -r requirements/requirements-tests.txt
-nodeenv==1.6.0
+netcdf4==1.6.2
+    # via
+    #   -r requirements/requirements-tests.txt
+    #   -r requirements/requirements.txt
+nodeenv==1.7.0
     # via pre-commit
-num2words==0.5.10
+num2words==0.5.12
     # via
     #   -r requirements/requirements.txt
     #   quantulum3
-numpy==1.22.3
+numpy==1.23.5
     # via
     #   -r requirements/requirements-tests.txt
     #   -r requirements/requirements.txt
     #   bokeh
+    #   cfgrib
     #   cftime
     #   dask
+    #   eccodes
     #   netcdf4
     #   pandas
     #   rasterio
     #   snuggs
     #   xarray
-openpyxl==3.0.9
+openpyxl==3.0.10
     # via
     #   -r requirements/requirements.txt
-    #   frictionless
+    #   hdx-python-utilities
 packaging==21.3
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements-tests.txt
     #   -r requirements/requirements.txt
     #   bokeh
+    #   build
     #   dask
     #   distributed
+    #   geopandas
     #   pytest
     #   rioxarray
     #   sphinx
     #   tox
     #   xarray
-pandas==1.4.2
+pandas==1.5.2
     # via
     #   -r requirements/requirements.txt
     #   dask
     #   geopandas
     #   xarray
-partd==1.2.0
+partd==1.3.0
     # via
     #   -r requirements/requirements.txt
     #   dask
-pbr==5.8.1
+pbr==5.11.0
     # via
     #   -r requirements/requirements-docs.txt
     #   sphinxcontrib-apidoc
-pep517==0.12.0
-    # via pip-tools
-petl==1.7.8
+pep517==0.13.0
+    # via build
+petl==1.7.12
     # via
     #   -r requirements/requirements.txt
     #   frictionless
-pillow==9.1.0
+pillow==9.3.0
     # via
     #   -r requirements/requirements.txt
     #   bokeh
-pip-tools==6.6.0
+pip-tools==6.10.0
     # via -r requirements/dev.in
-platformdirs==2.5.2
+platformdirs==2.5.4
     # via
     #   -r requirements/requirements-tests.txt
     #   virtualenv
@@ -368,16 +385,15 @@ pockets==0.9.1
     # via
     #   -r requirements/requirements.txt
     #   sphinxcontrib-napoleon
-pre-commit==2.18.1
+pre-commit==2.20.0
     # via -r requirements/dev.in
-psutil==5.9.0
+psutil==5.9.4
     # via
     #   -r requirements/requirements.txt
     #   distributed
 py==1.11.0
     # via
     #   -r requirements/requirements-tests.txt
-    #   pytest
     #   tox
 pyasn1==0.4.8
     # via
@@ -388,18 +404,22 @@ pycparser==2.21
     # via
     #   -r requirements/requirements.txt
     #   cffi
-pydantic==1.9.0
-    # via -r requirements/requirements.txt
-pygments==2.11.2
+pydantic==1.10.2
+    # via
+    #   -r requirements/requirements.txt
+    #   inflect
+pygments==2.13.0
     # via
     #   -r requirements/requirements-docs.txt
+    #   -r requirements/requirements.txt
+    #   rich
     #   sphinx
-pyopenssl==22.0.0
+pyopenssl==22.1.0
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-api
     #   ndg-httpsclient
-pyparsing==3.0.8
+pyparsing==3.0.9
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements-tests.txt
@@ -410,23 +430,23 @@ pyphonetics==0.5.3
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-country
-pyproj==3.3.0
+pyproj==3.4.0
     # via
     #   -r requirements/requirements.txt
     #   geopandas
     #   rioxarray
-pyrsistent==0.18.1
+pyrsistent==0.19.2
     # via
     #   -r requirements/requirements.txt
     #   jsonschema
-pytest==7.1.1
+pytest==7.2.0
     # via
     #   -r requirements/requirements-tests.txt
     #   pytest-cov
     #   pytest-mock
-pytest-cov==3.0.0
+pytest-cov==4.0.0
     # via -r requirements/requirements-tests.txt
-pytest-mock==3.7.0
+pytest-mock==3.10.0
     # via -r requirements/requirements-tests.txt
 python-dateutil==2.8.2
     # via
@@ -435,16 +455,16 @@ python-dateutil==2.8.2
     #   hdx-python-utilities
     #   libhxl
     #   pandas
-python-io-wrapper==0.2
+python-io-wrapper==0.3.1
     # via
     #   -r requirements/requirements.txt
     #   libhxl
-python-slugify==6.1.1
+python-slugify==7.0.0
     # via
     #   -r requirements/requirements.txt
     #   ckanapi
     #   frictionless
-pytz==2022.1
+pytz==2022.6
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements.txt
@@ -459,11 +479,11 @@ pyyaml==6.0
     #   frictionless
     #   pre-commit
     #   tableschema-to-template
-quantulum3==0.7.10
+quantulum3==0.7.11
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-api
-rasterio==1.2.10
+rasterio==1.3.4
     # via
     #   -r requirements/requirements.txt
     #   rioxarray
@@ -471,18 +491,18 @@ ratelimit==2.2.1
     # via
     #   -r requirements/requirements.txt
     #   hdx-python-utilities
-requests==2.27.1
+requests==2.28.1
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements.txt
+    #   cdsapi
     #   ckanapi
-    #   exchangerates
     #   frictionless
     #   libhxl
     #   requests-cache
     #   requests-file
     #   sphinx
-requests-cache==0.9.3
+requests-cache==0.9.7
     # via
     #   -r requirements/requirements.txt
     #   libhxl
@@ -494,21 +514,21 @@ rfc3986==2.0.0
     # via
     #   -r requirements/requirements.txt
     #   frictionless
-rioxarray==0.10.3
-    # via -r requirements/requirements.txt
-ruamel.yaml==0.17.21
+rich==12.6.0
     # via
     #   -r requirements/requirements.txt
-    #   hdx-python-utilities
-ruamel.yaml.clib==0.2.6
+    #   typer
+rioxarray==0.13.1
+    # via -r requirements/requirements.txt
+ruamel-yaml==0.17.21
     # via
     #   -r requirements/requirements.txt
-    #   ruamel.yaml
-shapely==1.8.1.post1
+    #   hdx-python-utilities
+shapely==1.8.5.post1
     # via
     #   -r requirements/requirements.txt
     #   geopandas
-shellingham==1.4.0
+shellingham==1.5.0
     # via
     #   -r requirements/requirements.txt
     #   typer
@@ -520,9 +540,7 @@ six==1.16.0
     # via
     #   -r requirements/requirements-tests.txt
     #   -r requirements/requirements.txt
-    #   basicauth
     #   ckanapi
-    #   exchangerates
     #   fiona
     #   isodate
     #   jsonpath-ng
@@ -533,8 +551,6 @@ six==1.16.0
     #   sphinxcontrib-napoleon
     #   tox
     #   url-normalize
-    #   validators
-    #   virtualenv
 snowballstemmer==2.2.0
     # via
     #   -r requirements/requirements-docs.txt
@@ -547,12 +563,12 @@ sortedcontainers==2.4.0
     # via
     #   -r requirements/requirements.txt
     #   distributed
-sphinx==4.5.0
+sphinx==5.3.0
     # via
     #   -r requirements/requirements-docs.txt
     #   sphinx-rtd-theme
     #   sphinxcontrib-apidoc
-sphinx-rtd-theme==1.0.0
+sphinx-rtd-theme==1.1.1
     # via -r requirements/requirements-docs.txt
 sphinxcontrib-apidoc==0.3.0
     # via -r requirements/requirements-docs.txt
@@ -591,8 +607,8 @@ stringcase==1.2.0
 tableschema-to-template==0.0.12
     # via
     #   -r requirements/requirements.txt
-    #   frictionless
-tabulate==0.8.9
+    #   hdx-python-utilities
+tabulate==0.9.0
     # via
     #   -r requirements/requirements.txt
     #   frictionless
@@ -605,17 +621,8 @@ text-unidecode==1.3
     #   -r requirements/requirements.txt
     #   python-slugify
 toml==0.10.2
-    # via
-    #   -r requirements/requirements-tests.txt
-    #   pre-commit
-    #   tox
-tomli==2.0.1
-    # via
-    #   -r requirements/requirements-tests.txt
-    #   coverage
-    #   pep517
-    #   pytest
-toolz==0.11.2
+    # via pre-commit
+toolz==0.12.0
     # via
     #   -r requirements/requirements.txt
     #   dask
@@ -626,23 +633,22 @@ tornado==6.1
     #   -r requirements/requirements.txt
     #   bokeh
     #   distributed
-tox==3.25.0
+tox==3.27.1
     # via -r requirements/requirements-tests.txt
-typer[all]==0.4.1
+tqdm==4.64.1
+    # via
+    #   -r requirements/requirements.txt
+    #   cdsapi
+typer[all]==0.7.0
     # via
     #   -r requirements/requirements.txt
     #   frictionless
-typing-extensions==4.2.0
+typing-extensions==4.4.0
     # via
     #   -r requirements/requirements.txt
     #   bokeh
     #   pydantic
-    #   typing-inspect
-typing-inspect==0.7.1
-    # via
-    #   -r requirements/requirements.txt
-    #   defopt
-unidecode==1.3.4
+unidecode==1.3.6
     # via
     #   -r requirements/requirements.txt
     #   libhxl
@@ -651,37 +657,37 @@ url-normalize==1.4.3
     # via
     #   -r requirements/requirements.txt
     #   requests-cache
-urllib3==1.26.9
+urllib3==1.26.13
     # via
     #   -r requirements/requirements-docs.txt
     #   -r requirements/requirements.txt
     #   distributed
     #   requests
     #   requests-cache
-validators==0.18.2
+validators==0.20.0
     # via
     #   -r requirements/requirements.txt
     #   frictionless
-virtualenv==20.14.1
+virtualenv==20.16.7
     # via
     #   -r requirements/requirements-tests.txt
     #   pre-commit
     #   tox
-wheel==0.37.1
+wheel==0.38.4
     # via
     #   -r requirements/requirements.txt
     #   libhxl
     #   pip-tools
-wrapt==1.14.0
+wrapt==1.14.1
     # via -r requirements/requirements.txt
-xarray[parallel]==2022.3.0
+xarray[parallel]==2022.11.0
     # via
     #   -r requirements/requirements.txt
     #   rioxarray
 xlrd==2.0.1
     # via
     #   -r requirements/requirements.txt
-    #   frictionless
+    #   hdx-python-utilities
 xlrd3==1.1.0
     # via
     #   -r requirements/requirements.txt
@@ -693,17 +699,11 @@ xlsxwriter==3.0.3
 xlwt==1.3.0
     # via
     #   -r requirements/requirements.txt
-    #   frictionless
-zict==2.1.0
+    #   hdx-python-utilities
+zict==2.2.0
     # via
     #   -r requirements/requirements.txt
     #   distributed
-zipp==3.8.0
-    # via
-    #   -r requirements/requirements-docs.txt
-    #   -r requirements/requirements.txt
-    #   importlib-metadata
-    #   importlib-resources
 
 # The following packages are considered to be unsafe in a requirements file:
 # pip
diff --git a/requirements/requirements-docs.txt b/requirements/requirements-docs.txt
index 654c40b2..9630d661 100644
--- a/requirements/requirements-docs.txt
+++ b/requirements/requirements-docs.txt
@@ -1,45 +1,49 @@
+#
+# This file is autogenerated by pip-compile with python 3.11
+# To update, run:
+#
+#    pip-compile --output-file=requirements/requirements-docs.txt requirements/docs.in
+#
 alabaster==0.7.12
     # via sphinx
-babel==2.10.1
+babel==2.11.0
     # via sphinx
-certifi==2021.10.8
+certifi==2022.9.24
     # via requests
-charset-normalizer==2.0.12
+charset-normalizer==2.1.1
     # via requests
 docutils==0.17.1
     # via
     #   sphinx
     #   sphinx-rtd-theme
-idna==3.3
+idna==3.4
     # via requests
-imagesize==1.3.0
+imagesize==1.4.1
     # via sphinx
-importlib-metadata==4.11.3
-    # via sphinx
-jinja2==3.1.1
+jinja2==3.1.2
     # via sphinx
 markupsafe==2.1.1
     # via jinja2
 packaging==21.3
     # via sphinx
-pbr==5.8.1
+pbr==5.11.0
     # via sphinxcontrib-apidoc
-pygments==2.11.2
+pygments==2.13.0
     # via sphinx
-pyparsing==3.0.8
+pyparsing==3.0.9
     # via packaging
-pytz==2022.1
+pytz==2022.6
     # via babel
-requests==2.27.1
+requests==2.28.1
     # via sphinx
 snowballstemmer==2.2.0
     # via sphinx
-sphinx==4.5.0
+sphinx==5.3.0
     # via
     #   -r requirements/docs.in
     #   sphinx-rtd-theme
     #   sphinxcontrib-apidoc
-sphinx-rtd-theme==1.0.0
+sphinx-rtd-theme==1.1.1
     # via -r requirements/docs.in
 sphinxcontrib-apidoc==0.3.0
     # via -r requirements/docs.in
@@ -55,7 +59,5 @@ sphinxcontrib-qthelp==1.0.3
     # via sphinx
 sphinxcontrib-serializinghtml==1.1.5
     # via sphinx
-urllib3==1.26.9
+urllib3==1.26.13
     # via requests
-zipp==3.8.0
-    # via importlib-metadata
diff --git a/requirements/requirements.txt b/requirements/requirements.txt
index ee21cd75..f522cd30 100644
--- a/requirements/requirements.txt
+++ b/requirements/requirements.txt
@@ -1,39 +1,57 @@
+#
+# This file is autogenerated by pip-compile with python 3.11
+# To update, run:
+#
+#    pip-compile --output-file=requirements/requirements.txt setup.cfg
+#
 affine==2.3.1
     # via rasterio
 appdirs==1.4.4
     # via requests-cache
-attrs==21.4.0
+attrs==22.1.0
     # via
     #   cattrs
+    #   cfgrib
+    #   eccodes
     #   fiona
     #   jsonlines
     #   jsonschema
     #   rasterio
     #   requests-cache
-basicauth==0.4.1
-    # via hdx-python-utilities
-bokeh==2.4.2
+bokeh==2.4.3
     # via dask
-cattrs==1.10.0
+cattrs==22.2.0
     # via requests-cache
-certifi==2021.10.8
+cdsapi==0.5.1
+    # via aa-toolbox (setup.cfg)
+certifi==2022.9.24
     # via
     #   fiona
     #   pyproj
     #   rasterio
     #   requests
-cffi==1.15.0
-    # via cryptography
-chardet==4.0.0
+cffi==1.15.1
+    # via
+    #   cryptography
+    #   eccodes
+cfgrib==0.9.10.3
+    # via aa-toolbox (setup.cfg)
+cftime==1.6.2
+    # via
+    #   aa-toolbox (setup.cfg)
+    #   netcdf4
+chardet==5.0.0
     # via frictionless
-charset-normalizer==2.0.12
+charset-normalizer==2.1.1
     # via requests
 ckanapi==4.7
     # via hdx-python-api
-click==8.1.2
+click==8.1.3
     # via
+    #   cfgrib
     #   click-plugins
     #   cligj
+    #   dask
     #   distributed
     #   fiona
     #   rasterio
@@ -46,15 +64,17 @@ cligj==0.7.2
     # via
     #   fiona
     #   rasterio
-cloudpickle==2.0.0
+cloudpickle==2.2.0
     # via
     #   dask
     #   distributed
-colorama==0.4.4
+colorama==0.4.6
     # via typer
-cryptography==36.0.2
+commonmark==0.9.1
+    # via rich
+cryptography==38.0.3
     # via pyopenssl
-dask==2022.04.1
+dask[complete]==2022.11.1
     # via
     #   distributed
     #   xarray
@@ -62,9 +82,9 @@ decorator==5.1.1
     # via
     #   jsonpath-ng
     #   validators
-defopt==6.3.0
+defopt==6.4.0
     # via hdx-python-api
-distributed==2022.04.1
+distributed==2022.11.1
     # via dask
 dnspython==2.2.1
     # via email-validator
@@ -74,106 +94,113 @@ docopt==0.6.2
     #   num2words
 docutils==0.17.1
     # via defopt
-email-validator==1.1.3
+eccodes==1.5.0
+    # via cfgrib
+email-validator==1.3.0
     # via hdx-python-api
 et-xmlfile==1.1.0
     # via openpyxl
-exchangerates==0.3.4
-    # via hdx-python-country
-fiona==1.8.21
+findlibs==0.0.2
+    # via eccodes
+fiona==1.8.22
     # via geopandas
-frictionless[excel,json]==4.34.0
+frictionless==4.40.8
     # via hdx-python-utilities
-fsspec==2022.3.0
+fsspec==2022.11.0
     # via dask
-geopandas==0.10.2
+geopandas==0.12.1
     # via aa-toolbox (setup.cfg)
-hdx-python-api==5.6.5
+hdx-python-api==5.9.2
     # via aa-toolbox (setup.cfg)
-hdx-python-country==3.2.1
+hdx-python-country==3.4.3
     # via
     #   aa-toolbox (setup.cfg)
     #   hdx-python-api
-hdx-python-utilities==3.3.1
+hdx-python-utilities==3.4.8
     # via hdx-python-country
 heapdict==1.0.1
     # via zict
-idna==3.3
+idna==3.4
     # via
     #   email-validator
     #   requests
 ijson==3.1.4
-    # via frictionless
-importlib-resources==5.7.1
-    # via jsonschema
-inflect==5.5.2
+    # via hdx-python-utilities
+inflect==6.0.2
     # via quantulum3
 isodate==0.6.1
     # via frictionless
-jinja2==3.1.1
+jinja2==3.1.2
     # via
     #   bokeh
     #   dask
     #   distributed
     #   frictionless
-jsonlines==3.0.0
-    # via frictionless
+jsonlines==3.1.0
+    # via hdx-python-utilities
 jsonpath-ng==1.5.3
     # via libhxl
-jsonschema==4.4.0
+jsonschema==4.17.1
     # via
     #   frictionless
     #   tableschema-to-template
-libhxl==4.25
+libhxl==4.27
     # via hdx-python-country
 locket==1.0.0
-    # via partd
+    # via
+    #   distributed
+    #   partd
 loguru==0.6.0
     # via hdx-python-utilities
-lxml==4.8.0
-    # via exchangerates
-marko==1.2.0
+makefun==1.15.0
+    # via hdx-python-api
+marko==1.2.2
     # via frictionless
 markupsafe==2.1.1
     # via jinja2
-msgpack==1.0.3
+msgpack==1.0.4
     # via distributed
 munch==2.5.0
     # via fiona
-mypy-extensions==0.4.3
-    # via typing-inspect
 ndg-httpsclient==0.5.1
     # via hdx-python-api
-num2words==0.5.10
+netcdf4==1.6.2
+    # via aa-toolbox (setup.cfg)
+num2words==0.5.12
     # via quantulum3
-numpy==1.22.3
+numpy==1.23.5
     # via
     #   aa-toolbox (setup.cfg)
     #   bokeh
+    #   cfgrib
+    #   cftime
     #   dask
+    #   eccodes
+    #   netcdf4
     #   pandas
     #   rasterio
     #   snuggs
     #   xarray
-openpyxl==3.0.9
-    # via frictionless
+openpyxl==3.0.10
+    # via hdx-python-utilities
 packaging==21.3
     # via
     #   bokeh
     #   dask
     #   distributed
+    #   geopandas
     #   rioxarray
     #   xarray
-pandas==1.4.2
+pandas==1.5.2
     # via
     #   dask
     #   geopandas
     #   xarray
-partd==1.2.0
+partd==1.3.0
     # via dask
-petl==1.7.8
+petl==1.7.12
     # via frictionless
-pillow==9.1.0
+pillow==9.3.0
     # via bokeh
 ply==3.11
     # via
@@ -181,7 +208,7 @@ ply==3.11
     #   libhxl
 pockets==0.9.1
     # via sphinxcontrib-napoleon
-psutil==5.9.0
+psutil==5.9.4
     # via distributed
 pyasn1==0.4.8
     # via
@@ -189,37 +216,42 @@ pyasn1==0.4.8
     #   ndg-httpsclient
 pycparser==2.21
     # via cffi
-pydantic==1.9.0
-    # via aa-toolbox (setup.cfg)
-pyopenssl==22.0.0
+pydantic==1.10.2
+    # via
+    #   aa-toolbox (setup.cfg)
+    #   inflect
+pygments==2.13.0
+    # via rich
+pyopenssl==22.1.0
     # via
     #   hdx-python-api
     #   ndg-httpsclient
-pyparsing==3.0.8
+pyparsing==3.0.9
     # via
     #   packaging
     #   snuggs
 pyphonetics==0.5.3
     # via hdx-python-country
-pyproj==3.3.0
+pyproj==3.4.0
     # via
     #   geopandas
     #   rioxarray
-pyrsistent==0.18.1
+pyrsistent==0.19.2
     # via jsonschema
 python-dateutil==2.8.2
     # via
+    #   aa-toolbox (setup.cfg)
     #   frictionless
     #   hdx-python-utilities
     #   libhxl
     #   pandas
-python-io-wrapper==0.2
+python-io-wrapper==0.3.1
     # via libhxl
-python-slugify==6.1.1
+python-slugify==7.0.0
     # via
     #   ckanapi
     #   frictionless
-pytz==2022.1
+pytz==2022.6
     # via pandas
 pyyaml==6.0
     # via
@@ -229,44 +261,44 @@ pyyaml==6.0
     #   distributed
     #   frictionless
     #   tableschema-to-template
-quantulum3==0.7.10
+quantulum3==0.7.11
     # via hdx-python-api
-rasterio==1.2.10
-    # via rioxarray
+rasterio==1.3.4
+    # via
+    #   aa-toolbox (setup.cfg)
+    #   rioxarray
 ratelimit==2.2.1
     # via hdx-python-utilities
-requests==2.27.1
+requests==2.28.1
     # via
     #   aa-toolbox (setup.cfg)
+    #   cdsapi
     #   ckanapi
-    #   exchangerates
     #   frictionless
     #   libhxl
     #   requests-cache
     #   requests-file
-requests-cache==0.9.3
+requests-cache==0.9.7
     # via libhxl
 requests-file==1.5.1
     # via hdx-python-utilities
 rfc3986==2.0.0
     # via frictionless
-rioxarray==0.10.3
+rich==12.6.0
+    # via typer
+rioxarray==0.13.1
     # via aa-toolbox (setup.cfg)
-ruamel.yaml==0.17.21
+ruamel-yaml==0.17.21
     # via hdx-python-utilities
-ruamel.yaml.clib==0.2.6
-    # via ruamel.yaml
-shapely==1.8.1.post1
+shapely==1.8.5.post1
     # via geopandas
-shellingham==1.4.0
+shellingham==1.5.0
     # via typer
 simpleeval==0.9.12
     # via frictionless
 six==1.16.0
     # via
-    #   basicauth
     #   ckanapi
-    #   exchangerates
     #   fiona
     #   isodate
     #   jsonpath-ng
@@ -276,7 +308,6 @@ six==1.16.0
     #   requests-file
     #   sphinxcontrib-napoleon
     #   url-normalize
-    #   validators
 snuggs==1.4.7
     # via rasterio
 sortedcontainers==2.4.0
@@ -286,14 +317,14 @@ sphinxcontrib-napoleon==0.7
 stringcase==1.2.0
     # via frictionless
 tableschema-to-template==0.0.12
-    # via frictionless
-tabulate==0.8.9
+    # via hdx-python-utilities
+tabulate==0.9.0
     # via frictionless
 tblib==1.7.0
     # via distributed
 text-unidecode==1.3
     # via python-slugify
-toolz==0.11.2
+toolz==0.12.0
     # via
     #   dask
     #   distributed
@@ -302,49 +333,45 @@ tornado==6.1
     # via
     #   bokeh
     #   distributed
-typer[all]==0.4.1
+tqdm==4.64.1
+    # via cdsapi
+typer[all]==0.7.0
     # via frictionless
-typing-extensions==4.2.0
+typing-extensions==4.4.0
     # via
-    #   aa-toolbox (setup.cfg)
     #   bokeh
     #   pydantic
-    #   typing-inspect
-typing-inspect==0.7.1
-    # via defopt
-unidecode==1.3.4
+unidecode==1.3.6
     # via
     #   libhxl
     #   pyphonetics
 url-normalize==1.4.3
     # via requests-cache
-urllib3==1.26.9
+urllib3==1.26.13
     # via
     #   distributed
     #   requests
     #   requests-cache
-validators==0.18.2
+validators==0.20.0
     # via frictionless
-wheel==0.37.1
+wheel==0.38.4
     # via libhxl
-wrapt==1.14.0
+wrapt==1.14.1
     # via aa-toolbox (setup.cfg)
-xarray[parallel]==2022.3.0
+xarray[parallel]==2022.11.0
     # via
     #   aa-toolbox (setup.cfg)
     #   rioxarray
 xlrd==2.0.1
-    # via frictionless
+    # via hdx-python-utilities
 xlrd3==1.1.0
     # via libhxl
 xlsxwriter==3.0.3
     # via tableschema-to-template
 xlwt==1.3.0
-    # via frictionless
-zict==2.1.0
+    # via hdx-python-utilities
+zict==2.2.0
     # via distributed
-zipp==3.8.0
-    # via importlib-resources
 
 # The following packages are considered to be unsafe in a requirements file:
 # setuptools
diff --git a/setup.cfg b/setup.cfg
index f512129a..f045140e 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -27,19 +27,27 @@ classifiers =
 packages = find:
 package_dir = = src
 include_package_data = true
-python_requires = >= 3.6
-# TODO: Remove the rasterio requirements once 1.3 is released
+python_requires = >= 3.7
+# TODO: split up this list into custom installs
+# TODO: Remove importlib requirement once error fixed
+#  https://stackoverflow.com/questions/73929564/
+#  entrypoints-object-has-no-attribute-get-digital-ocean
 install_requires =
+    cdsapi
+    cfgrib
+    cftime
     geopandas
     hdx-python-api>=5.6.4
     hdx-python-country
+    importlib-metadata==4.13.0;python_version=='3.7'
+    netCDF4
     numpy
     pydantic
+    python-dateutil
     pyyaml
+    rasterio
     requests
-    rasterio>=1.3a4;python_version=='3.10'
     rioxarray
-    typing-extensions
     wrapt
     xarray[parallel]
 
diff --git a/src/aatoolbox/__init__.py b/src/aatoolbox/__init__.py
index 5dd8f502..9f1721c9 100644
--- a/src/aatoolbox/__init__.py
+++ b/src/aatoolbox/__init__.py
@@ -5,26 +5,51 @@
 anticipatory action.
 """
 from aatoolbox.config.countryconfig import (
+    CountryConfig,
     create_country_config,
     create_custom_country_config,
 )
+from aatoolbox.datasources.chirps.chirps import ChirpsDaily, ChirpsMonthly
 from aatoolbox.datasources.codab.codab import CodAB
 from aatoolbox.datasources.fewsnet.fewsnet import FewsNet
+from aatoolbox.datasources.glofas.forecast import (
+    GlofasForecast,
+    GlofasReforecast,
+)
+from aatoolbox.datasources.glofas.reanalysis import GlofasReanalysis
 from aatoolbox.datasources.iri.iri_seasonal_forecast import (
     IriForecastDominant,
     IriForecastProb,
 )
+from aatoolbox.datasources.usgs.ndvi_products import (
+    UsgsNdviMedianAnomaly,
+    UsgsNdviPctMedian,
+    UsgsNdviSmoothed,
+    UsgsNdviYearDifference,
+)
+from aatoolbox.utils import raster
 from aatoolbox.utils.geoboundingbox import GeoBoundingBox
 
 from ._version import version as __version__
 
 __all__ = (
+    "CountryConfig",
     "create_country_config",
     "create_custom_country_config",
+    "ChirpsDaily",
+    "ChirpsMonthly",
     "CodAB",
-    "GeoBoundingBox",
+    "FewsNet",
+    "GlofasReanalysis",
+    "GlofasForecast",
+    "GlofasReforecast",
     "IriForecastProb",
     "IriForecastDominant",
-    "FewsNet",
+    "UsgsNdviSmoothed",
+    "UsgsNdviPctMedian",
+    "UsgsNdviMedianAnomaly",
+    "UsgsNdviYearDifference",
+    "GeoBoundingBox",
+    "raster",
     "__version__",
 )
diff --git a/src/aatoolbox/config/countries/bfa.yaml b/src/aatoolbox/config/countries/bfa.yaml
index 869b7dd5..ad9b975a 100644
--- a/src/aatoolbox/config/countries/bfa.yaml
+++ b/src/aatoolbox/config/countries/bfa.yaml
@@ -5,3 +5,37 @@ codab:
   admin_level_max: 3
 fewsnet:
   region_name: west-africa
+glofas:
+  reporting_points:
+  - name: BOUGOURI-BA a DIEBOUGOU
+    lon: -3.15
+    lat: 10.95
+  - name: Bassieri
+    lon: 0.35
+    lat: 12.75
+  - name: Faga A Liptougou
+    lon: -0.45
+    lat: 13.25
+  - name: Liptougou
+    lon: 0.35
+    lat: 13.25
+  - name: Manni A Manni
+    lon: -0.35
+    lat: 13.15
+  - name: Sebba
+    lon: 0.45
+    lat: 13.45
+  - name: Sourou A Pont De Leri Nord
+    lon: -3.55
+    lat: 12.65
+  - name: VOLTA ROUGE a DAKAYE
+    lon: -1.65
+    lat: 11.75
+  - name: Volta Blanche A Wayen
+    lon: -1.15
+    lat: 12.35
+  - name: Volta Noire A Samandeni
+    lon: -4.45
+    lat: 11.45
+usgs_ndvi:
+  area_name: west-africa
diff --git a/src/aatoolbox/config/countries/bgd.yaml b/src/aatoolbox/config/countries/bgd.yaml
index 351a645d..e3794941 100644
--- a/src/aatoolbox/config/countries/bgd.yaml
+++ b/src/aatoolbox/config/countries/bgd.yaml
@@ -3,3 +3,35 @@ codab:
   hdx_dataset_name: bgd_adm_bbs_20201113_SHP.zip
   layer_base_name: bgd_adm_bbs_20201113_SHP/bgd_admbnda_adm{admin_level}_bbs_20201113.shp
   admin_level_max: 4
+glofas:
+  reporting_points:
+  - name: Bahadurabad
+    lon: 89.65
+    lat: 25.15
+  - name: Bhairab Bazar
+    lon: 90.95
+    lat: 24.05
+  - name: Hardinge Bridge
+    lon: 89.05
+    lat: 24.05
+  - name: Jariajanjail
+    lon: 90.65
+    lat: 24.95
+  - name: Kanaighat
+    lon: 92.15
+    lat: 25.05
+  - name: Kaunia
+    lon: 89.45
+    lat: 25.75
+  - name: Laurergarh Saktiarkhola
+    lon: 91.25
+    lat: 25.15
+  - name: Na
+    lon: 89.75
+    lat: 25.55
+  - name: Sherpur
+    lon: 91.65
+    lat: 24.65
+  - name: Sunamganj
+    lon: 91.45
+    lat: 25.05
diff --git a/src/aatoolbox/config/countries/cod.yaml b/src/aatoolbox/config/countries/cod.yaml
index 8d75b347..8dd9e9ce 100644
--- a/src/aatoolbox/config/countries/cod.yaml
+++ b/src/aatoolbox/config/countries/cod.yaml
@@ -3,8 +3,52 @@ codab:
   hdx_dataset_name: cod_admbnda_rgc_itos_20190911_SHP.zip
   layer_base_name: cod_admbnda_adm{admin_level}_rgc_itos_20190911.shp
   admin_level_max: 1
-  # this is the name of the adm2 layer
-  custom_layer_names:
-    - cod_admbnda_adm2_rgc_20190911.shp
+  admin2_name: cod_admbnda_adm2_rgc_20190911.shp
 fewsnet:
   region_name: southern-africa
+glofas:
+  reporting_points:
+  - name: Bangassou
+    lon: 22.85
+    lat: 4.75
+  - name: Congo A Kinshasa-Est-H-/Inga-Q-
+    lon: 15.25
+    lat: -4.35
+  - name: Congo River At Kindu
+    lon: 25.95
+    lat: -2.95
+  - name: Kasai River At Ilebo -Ex- Port Francqui-
+    lon: 20.55
+    lat: -4.35
+  - name: Kinshasa
+    lon: 15.25
+    lat: -4.35
+  - name: Lualaba River At Kisangani -Ex- Stanleyville-
+    lon: 25.15
+    lat: 0.55
+  - name: Lufubu At Chipili
+    lon: 29.05
+    lat: -10.75
+  - name: Luongo At Mwendakashiba
+    lon: 29.05
+    lat: -10.45
+  - name: Manovo Manovo
+    lon: 22.95
+    lat: 8.15
+  - name: Mbomou Bangassou
+    lon: 22.85
+    lat: 4.75
+  - name: Na
+    lon: 27.45
+    lat: -10.45
+  - name: Oubangui Mobaye
+    lon: 21.25
+    lat: 4.35
+  - name: Semliki 85205
+    lon: 29.85
+    lat: 0.75
+  - name: Zemio
+    lon: 25.15
+    lat: 5.05
+usgs_ndvi:
+  area_name: southern-africa
diff --git a/src/aatoolbox/config/countries/eth.yaml b/src/aatoolbox/config/countries/eth.yaml
index 80512e57..4b530d25 100644
--- a/src/aatoolbox/config/countries/eth.yaml
+++ b/src/aatoolbox/config/countries/eth.yaml
@@ -2,6 +2,86 @@ iso3: eth
 codab:
   hdx_dataset_name: ETH Administrative Divisions Shapefiles.zip
   layer_base_name: eth_admbnda_adm{admin_level}_csa_bofed_20201008.shp
-  admin_level_max: 2
+  admin_level_max: 3
+  admin0_name: eth_admbnda_adm0_csa_bofedb_itos_2021
 fewsnet:
   region_name: east-africa
+glofas:
+  reporting_points:
+  - name: Abetlti
+    lon: 37.55
+    lat: 8.25
+  - name: Adaitu
+    lon: 40.65
+    lat: 11.15
+  - name: Asetndabo
+    lon: 37.15
+    lat: 7.75
+  - name: Awash Wenz At Hombole
+    lon: 38.85
+    lat: 8.35
+  - name: Awash Wenz At Melka Kuntire
+    lon: 38.55
+    lat: 8.75
+  - name: Bahir Dar
+    lon: 37.35
+    lat: 11.55
+  - name: Eddiem
+    lon: 34.95
+    lat: 11.25
+  - name: Gambella
+    lon: 34.65
+    lat: 8.15
+  - name: Melka Sedi
+    lon: 40.15
+    lat: 9.15
+  - name: Metahara
+    lon: 39.95
+    lat: 8.85
+  - name: Na 1
+    lon: 42.05
+    lat: 4.75
+  - name: Na 2
+    lon: 44.85
+    lat: 5.05
+  - name: Na 3
+    lon: 44.45
+    lat: 5.45
+  - name: Na 4
+    lon: 44.05
+    lat: 5.65
+  - name: Na 5
+    lon: 43.05
+    lat: 6.05
+  - name: Nr. Addis Zemen
+    lon: 37.75
+    lat: 12.05
+  - name: Nr. Bahir Dar
+    lon: 37.65
+    lat: 11.85
+  - name: Nr. Kessie
+    lon: 38.15
+    lat: 10.05
+  - name: Nr. Marawi
+    lon: 37.05
+    lat: 11.35
+  - name: Omo River Delta
+    lon: 36.05
+    lat: 4.85
+  - name: Sekta
+    lon: 36.75
+    lat: 7.65
+  - name: Shetbe
+    lon: 36.45
+    lat: 7.45
+  - name: Tendaho
+    lon: 40.95
+    lat: 11.65
+  - name: Woltkite 1
+    lon: 37.75
+    lat: 8.25
+  - name: Woltkite 2
+    lon: 37.55
+    lat: 8.35
+usgs_ndvi:
+  area_name: east-africa
diff --git a/src/aatoolbox/config/countries/mwi.yaml b/src/aatoolbox/config/countries/mwi.yaml
index fcc1641e..e9549d59 100644
--- a/src/aatoolbox/config/countries/mwi.yaml
+++ b/src/aatoolbox/config/countries/mwi.yaml
@@ -5,3 +5,43 @@ codab:
   admin_level_max: 3
 fewsnet:
   region_name: southern-africa
+glofas:
+  reporting_points:
+  - name: Chikwawa
+    lon: 34.85
+    lat: -16.05
+  - name: Lichenya River At Milonde -65312108-
+    lon: 35.45
+    lat: -16.05
+  - name: Liwonde
+    lon: 35.25
+    lat: -15.05
+  - name: Luchenza River At Luchenza -65312110-
+    lon: 35.25
+    lat: -15.95
+  - name: Luwewya At Majikapotwe/Zayuka -65312502-
+    lon: 34.15
+    lat: -11.75
+  - name: Mzimba At Muweru Bukuluti -65312404-
+    lon: 33.65
+    lat: -11.95
+  - name: Na 1
+    lon: 35.05
+    lat: -16.45
+  - name: Na 2
+    lon: 34.95
+    lat: -16.25
+  - name: Phewzi 65312414
+    lon: 33.95
+    lat: -10.95
+  - name: Rumphi At Chiweta -65312407-
+    lon: 34.05
+    lat: -10.65
+  - name: Runyina River At Mjuma -65312406-
+    lon: 33.75
+    lat: -10.85
+  - name: S53 Roadbridge 65312601
+    lon: 34.15
+    lat: -12.85
+usgs_ndvi:
+  area_name: southern-africa
diff --git a/src/aatoolbox/config/countries/npl.yaml b/src/aatoolbox/config/countries/npl.yaml
index 1ebcde1c..837d4e27 100644
--- a/src/aatoolbox/config/countries/npl.yaml
+++ b/src/aatoolbox/config/countries/npl.yaml
@@ -3,5 +3,89 @@ codab:
   hdx_dataset_name: npl_admbnda_nd_20201117_SHP.zip
   layer_base_name: npl_admbnda_adm{admin_level}_nd_20201117.shp
   admin_level_max: 2
+  # Nepal also has districts that don't fall into a particular admin level
   custom_layer_names:
     - npl_admbnda_districts_nd_20201117
+glofas:
+  reporting_points:
+  - name: After Lothar Confluence
+    lon: 84.65
+    lat: 27.55
+  - name: Angsing
+    lon: 83.75
+    lat: 27.85
+  - name: Arughat
+    lon: 84.85
+    lat: 28.05
+  - name: Asaraghat
+    lon: 81.45
+    lat: 28.95
+  - name: Betrawati
+    lon: 85.15
+    lat: 28.05
+  - name: Bimalnagar
+    lon: 84.45
+    lat: 28.05
+  - name: Butwal
+    lon: 83.45
+    lat: 27.75
+  - name: Chatara
+    lon: 87.15
+    lat: 26.85
+  - name: Chepang
+    lon: 81.75
+    lat: 28.35
+  - name: Chisapani
+    lon: 81.25
+    lat: 28.65
+  - name: Devghat
+    lon: 84.45
+    lat: 27.75
+  - name: Dipayal
+    lon: 80.95
+    lat: 29.25
+  - name: Jomsom
+    lon: 83.65
+    lat: 28.75
+  - name: Kampughat
+    lon: 86.85
+    lat: 26.85
+  - name: Khurkot
+    lon: 86.05
+    lat: 27.25
+  - name: Kusum
+    lon: 82.15
+    lat: 27.95
+  - name: Mainachuli
+    lon: 87.85
+    lat: 26.65
+  - name: Majhitar 1
+    lon: 87.65
+    lat: 27.25
+  - name: Majhitar 2
+    lon: 87.45
+    lat: 26.95
+  - name: Pachuwarghat
+    lon: 85.75
+    lat: 27.55
+  - name: Parigaun
+    lon: 80.25
+    lat: 29.15
+  - name: Rabuwabazar
+    lon: 86.65
+    lat: 27.25
+  - name: Rai Goan
+    lon: 85.45
+    lat: 27.25
+  - name: Samajhighat
+    lon: 81.65
+    lat: 28.55
+  - name: Shisaghat
+    lon: 84.25
+    lat: 28.15
+  - name: Simle
+    lon: 87.15
+    lat: 27.05
+  - name: Turkighat
+    lon: 87.15
+    lat: 27.35
diff --git a/src/aatoolbox/config/countryconfig.py b/src/aatoolbox/config/countryconfig.py
index 13d287ac..59137a0c 100644
--- a/src/aatoolbox/config/countryconfig.py
+++ b/src/aatoolbox/config/countryconfig.py
@@ -1,6 +1,6 @@
 """Country configuration setting base class."""
 from pathlib import Path
-from typing import Dict, Optional, Union
+from typing import Dict, List, Optional, Tuple, Union
 
 from pydantic import BaseModel, root_validator, validator
 
@@ -20,7 +20,11 @@ class CodABConfig(BaseModel):
         change by a single custom_layer_number depending on the level. Should
         contain {admin_level} in place of the custom_layer_number.
     admin_level_max: int
-        The maximum admin level available in the shapefile.
+         The maximum admin level available in the shapefile, cannot be
+         greater than 4
+    admin{level}_name: str, optional
+        The names of any admin level layers that do not conform to the
+        layer_base_name pattern, where {level} ranges from 0 to 4
     custom_layer_names: list, optional
         Any additional layer names that don't fit into the admin level paradigm
     """
@@ -28,17 +32,43 @@ class CodABConfig(BaseModel):
     hdx_dataset_name: str
     layer_base_name: str
     admin_level_max: int
+    admin0_name: Optional[str]
+    admin1_name: Optional[str]
+    admin2_name: Optional[str]
+    admin3_name: Optional[str]
+    admin4_name: Optional[str]
     custom_layer_names: Optional[list]
 
     @validator("layer_base_name")
     def _validate_layer_base_name(cls, layer_base_name):
+        """Ensure that the layer basename contains {admin_level}."""
         if "{admin_level}" not in layer_base_name:
             raise ValueError(
-                "The layer base name must contain "
-                "an {admin_level} placeholder."
+                "In the COD AB section of the country configuration file, "
+                "layer_base_name must contain an {admin_level} placeholder."
             )
         return layer_base_name
 
+    @validator("admin_level_max")
+    def _validate_admin_level_max(cls, admin_level_max):
+        """Ensure that admin_level_max is between 0 and 4."""
+        if not 0 <= admin_level_max <= 4:
+            raise ValueError(
+                "In the COD AB section of the country configuration file, "
+                "admin_level_max must be between 0 and 4."
+            )
+        return admin_level_max
+
+    @root_validator(pre=False, skip_on_failure=True)
+    def _set_admin_levels(cls, values) -> dict:
+        """Set admin levels names using layer base name."""
+        for admin_level in range(values["admin_level_max"] + 1):
+            if values.get(f"admin{admin_level}_name", None) is None:
+                values[f"admin{admin_level}_name"] = values[
+                    "layer_base_name"
+                ].format(admin_level=admin_level)
+        return values
+
 
 class FewsNetConfig(BaseModel):
     """FEWS NET configuration.
@@ -81,15 +111,88 @@ def _set_region_code(cls, values) -> dict:
         return values
 
 
+class ReportingPoints(BaseModel):
+    """Coordinates of GloFAS reporting points."""
+
+    name: str
+    lon: float
+    lat: float
+
+
+class GlofasConfig(BaseModel):
+    """GloFAS configuration."""
+
+    reporting_points: List[ReportingPoints]
+
+
+class UsgsNdviConfig(BaseModel):
+    """USGS NDVI configuration.
+
+    Parameters
+    ----------
+    area_name : str
+        Name of the USGS NDVI coverage area the country belongs to.
+        Needed to download the regional NDVI data.
+    """
+
+    area_name_mapping: Dict[str, Tuple[str, str]] = {
+        "north-africa": ("africa/north", "na"),
+        "east-africa": ("africa/east", "ea"),
+        "southern-africa": ("africa/southern", "sa"),
+        "west-africa": ("africa/west", "wa"),
+        "central-asia": ("asia/centralasia", "cta"),
+        "yemen": ("asia/middleeast/yemen", "yem"),
+        "central-america": ("lac/camcar/centralamerica", "ca"),
+        "hispaniola": ("lac/camcar/caribbean/hispaniola", "hi"),
+    }
+    area_name: str
+
+    @validator("area_name")
+    def area_name_valid(cls, v, values) -> str:
+        """Check that area_name is valid."""
+        valid_area_names = values["area_name_mapping"].keys()
+        if v not in valid_area_names:
+            raise ValueError(
+                f"Invalid area name: {v}. "
+                f"Should be one of {', '.join(valid_area_names)}"
+            )
+        return v
+
+    @root_validator(pre=False, skip_on_failure=True)
+    def _set_area_codes(cls, values) -> dict:
+        """Set NDVI url and prefix from area."""
+        values["area_url"], values["area_prefix"] = values[
+            "area_name_mapping"
+        ][values["area_name"]]
+        return values
+
+
 class CountryConfig(BaseModel):
-    """Country configuration."""
+    """Country configuration.
+
+    Parameters
+    ----------
+    iso3 : str
+        Country ISO3, must be exactly 3 letters long
+    codab: CodABConfig, optional
+        Configuration object for COD AB
+    fewsnet: FewsNetConfig, optional
+        Configuration object for FEWS NET
+    glofas: GlofasConfig, optional
+        Configuration object for GloFAS
+    usgs_ndvi: UsgsNdviConfig, optional
+        Configuration object for USGS NDVI
+    """
 
     iso3: str
     codab: Optional[CodABConfig]
     fewsnet: Optional[FewsNetConfig]
+    glofas: Optional[GlofasConfig]
+    usgs_ndvi: Optional[UsgsNdviConfig]
 
     @validator("iso3")
     def _validate_iso3(cls, iso3):
+        """Ensure ISO3 is length three and alphabet chars only."""
         return _validate_iso3(iso3)
 
 
@@ -138,6 +241,7 @@ def create_custom_country_config(filepath: Union[str, Path]) -> CountryConfig:
 
 
 def _validate_iso3(iso3: str):
+    """Ensure ISO3 is length three and alphabet chars only."""
     if len(iso3) != 3 or not str.isalpha(iso3):
         raise ValueError("ISO3 must be a three letter string.")
     return iso3.lower()
diff --git a/src/aatoolbox/datasources/chirps/__init__.py b/src/aatoolbox/datasources/chirps/__init__.py
new file mode 100644
index 00000000..51e9ea90
--- /dev/null
+++ b/src/aatoolbox/datasources/chirps/__init__.py
@@ -0,0 +1,1 @@
+"""Functionality to process CHIRPS data."""
diff --git a/src/aatoolbox/datasources/chirps/chirps.py b/src/aatoolbox/datasources/chirps/chirps.py
new file mode 100644
index 00000000..935e8c29
--- /dev/null
+++ b/src/aatoolbox/datasources/chirps/chirps.py
@@ -0,0 +1,623 @@
+"""Class to download and load CHIRPS observational precipitation data."""
+import calendar
+import logging
+import ssl
+from abc import abstractmethod
+from datetime import date, datetime
+from pathlib import Path
+from typing import Optional, Union
+from urllib.request import urlopen
+
+import cftime
+import pandas as pd
+import requests
+import xarray as xr
+
+from aatoolbox.config.countryconfig import CountryConfig
+from aatoolbox.datasources.datasource import DataSource
+from aatoolbox.utils.check_file_existence import check_file_existence
+from aatoolbox.utils.dates import get_date_from_user_input
+from aatoolbox.utils.geoboundingbox import GeoBoundingBox
+
+logger = logging.getLogger(__name__)
+
+_FIRST_AVAILABLE_DATE = date(year=1981, month=1, day=1)
+
+_VALID_RESOLUTIONS = (0.05, 0.25)
+
+_BASE_URL = "https://iridl.ldeo.columbia.edu/SOURCES/.UCSB/.CHIRPS/.v2p0/"
+
+
+class _Chirps(DataSource):
+    """
+    Base class object to retrieve CHIRPS observational precipitation data.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    geo_bounding_box: GeoBoundingBox
+        the bounding coordinates of the area that should be included in the
+        data.
+    frequency: str
+        time resolution of the data to be downloaded. It can be "daily" or
+        "monthly".
+    resolution: float, default = 0.05
+        resolution of data to be downloaded. Can be 0.05 or 0.25.
+    start_date: Optional[Union[datetime.date, str]], default = None
+        Data will be considered starting from date `start_date`.
+        Input can be an ISO8601 string or `datetime.date` object.
+        If None, it is set to 1981-1-1.
+    end_date: Optional[Union[datetime.date, str]], default = None
+        Data will be considered up to date `end_date`.
+        Input can be an ISO8601 string or `datetime.date` object.
+        If None, it is set to the date for which most recent data
+        is available.
+    """
+
+    @abstractmethod
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        geo_bounding_box: GeoBoundingBox,
+        frequency: str,
+        date_range_freq: str,
+        resolution: float = 0.05,
+        start_date: Optional[Union[date, str]] = None,
+        end_date: Optional[Union[date, str]] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            datasource_base_dir="chirps",
+            is_public=True,
+        )
+
+        # round coordinates to correspond with the grid the CHIRPS data is on
+        # non-rounded coordinates can be given to the URL which then
+        # automatically rounds them, but for file saving we prefer to do
+        # this ourselves
+        self._geobb = geo_bounding_box.round_coords(round_val=resolution)
+        self._frequency = frequency
+        self._resolution = resolution
+        self._date_range_freq = date_range_freq
+
+        if start_date is None:
+            start_date = _FIRST_AVAILABLE_DATE
+        if end_date is None:
+            end_date = self._get_last_available_date()
+
+        self._start_date = get_date_from_user_input(start_date)
+        self._end_date = get_date_from_user_input(end_date)
+
+        self._check_dates_validity()
+
+        if resolution not in _VALID_RESOLUTIONS:
+            raise ValueError(
+                f"The given resolution is {self._resolution}, which is "
+                "not available. Has to be 0.05 or 0.25."
+            )
+
+    def download(  # type: ignore
+        self,
+        clobber: bool = False,
+    ):
+        """
+        Download the CHIRPS observed precipitation as NetCDF file.
+
+        Parameters
+        ----------
+        clobber : bool, default = False
+            If True, overwrites existing raw files
+
+        Returns
+        -------
+        The folder where the data is downloaded.
+        """
+        # Create a list of date tuples
+        date_list = self._create_date_list(logging_level=logging.INFO)
+
+        # Data download
+        for d in date_list:
+            last_filepath = self._download_prep(d=d, clobber=clobber)
+
+        return last_filepath.parents[0]
+
+    def process(self, clobber: bool = False):
+        """
+        Process the CHIRPS data.
+
+        Should only be called after data has been download.
+
+        Parameters
+        ----------
+        clobber : bool, default = False
+            If True, overwrites existing processed files.
+
+        Returns
+        -------
+        The folder where the data is processed.
+        """
+        # Create a list with all raw data downloaded
+        filepath_list = self._get_to_be_processed_path_list()
+
+        for filepath in filepath_list:
+            try:
+                ds = xr.open_dataset(filepath, decode_times=False)
+            except ValueError as err:
+                raise ValueError(
+                    f"The dataset {filepath} is not a valid netcdf file: "
+                    "something probbly went wrong during the download. "
+                    "Try downloading the file again."
+                ) from err
+            processed_file_path = self._get_processed_path(filepath)
+            processed_file_path.parent.mkdir(parents=True, exist_ok=True)
+            last_filepath = self._process(
+                filepath=processed_file_path, ds=ds, clobber=clobber
+            )
+
+        return last_filepath.parents[0]
+
+    def load(self) -> xr.Dataset:
+        """
+        Load the CHIRPS data.
+
+        Should only be called after the data
+        has been downloaded and processed.
+
+        Returns
+        -------
+        The processed CHIRPS dataset.
+        """
+        # Get list of filepaths of files to be loaded
+
+        filepath_list = self._get_to_be_loaded_path_list()
+
+        # Merge all files in one dataset
+        if not filepath_list:
+            raise FileNotFoundError(
+                "Cannot find any netcdf file for the chosen combination "
+                "of frequency, resolution and area. Make sure "
+                "sure that you have already called the 'process' method."
+            )
+
+        try:
+            ds = xr.open_mfdataset(
+                filepath_list,
+            )
+            # include the names of all files that are included in the ds
+            ds.attrs["included_files"] = [f.stem for f in filepath_list]
+        except FileNotFoundError as err:
+            raise FileNotFoundError(
+                "Cannot find one or more netcdf files corresponding "
+                "to the selected range. Make sure that you already "
+                "downloaded and processed those data."
+            ) from err
+
+        return ds.rio.write_crs("EPSG:4326", inplace=True)
+
+    def _download_prep(self, d, clobber):
+        # Preparatory steps for the actual download
+        year = str(d.year)
+        month = f"{d.month:02d}"
+        day = f"{d.day:02d}"
+        output_filepath = self._get_raw_path(year=year, month=month, day=day)
+        output_filepath.parent.mkdir(parents=True, exist_ok=True)
+        url = self._get_url(year=year, month=month, day=day)
+        # Actual download
+        return self._download(
+            filepath=output_filepath,
+            url=url,
+            clobber=clobber,
+        )
+
+    def _create_date_list(self, logging_level=logging.DEBUG):
+        """Create list of tuples containing the range of dates of interest."""
+        date_list = pd.date_range(
+            self._start_date, self._end_date, freq=self._date_range_freq
+        ).tolist()
+
+        # Create a message containing information on the downloaded data
+        msg = (
+            f"{self._frequency.capitalize()} "
+            "data will be downloaded, starting from "
+            f"{date_list[0]} to {date_list[-1]}."
+        )
+
+        logging.log(logging_level, msg)
+
+        return date_list
+
+    def _check_dates_validity(self):
+        """Check dates vailidity."""
+        end_avail_date = self._get_last_available_date()
+
+        if (
+            not _FIRST_AVAILABLE_DATE
+            <= self._start_date
+            <= self._end_date
+            <= end_avail_date
+        ):
+            raise ValueError(
+                "Make sure that the input dates are ordered in the following "
+                f"way: {_FIRST_AVAILABLE_DATE} <= {self._start_date} <= "
+                f"{self._end_date} <= {end_avail_date}. The two dates above "
+                "indicate the range for which CHIRPS data are "
+                "currently available."
+            )
+
+    def _get_file_name_base(
+        self,
+        year: str,
+        month: str,
+    ) -> str:
+        if len(month) == 1:
+            month = f"0{month}"
+
+        file_name_base = (
+            f"{self._country_config.iso3}_chirps_"
+            f"{self._frequency}_{year}_{month}_"
+        )
+
+        return file_name_base
+
+    @abstractmethod
+    def _get_file_name(
+        self,
+        year: str,
+        month: str,
+        day: str,
+    ) -> str:
+        pass
+
+    def _get_raw_path(self, year: str, month: str, day: str) -> Path:
+        return self._raw_base_dir / self._get_file_name(
+            year=year, month=month, day=day
+        )
+
+    def _get_processed_path(self, raw_path: Path) -> Path:
+        return self._processed_base_dir / raw_path.parts[-1]
+
+    def _get_location_url(self):
+
+        location_url = (
+            f"X/%28{self._geobb.lon_min}%29%28{self._geobb.lon_max}"
+            f"%29RANGEEDGES/"
+            f"Y/%28{self._geobb.lat_max}%29%28{self._geobb.lat_min}"
+            f"%29RANGEEDGES/"
+        )
+
+        return location_url
+
+    @abstractmethod
+    def _get_url(self, year: str, month: str, day: str):
+        pass
+
+    @abstractmethod
+    def _get_last_available_date(self):
+        pass
+
+    def _get_to_be_processed_path_list(self):
+        """Get list of filepaths of files to be processed."""
+        date_list = self._create_date_list()
+
+        filepath_list = [
+            self._get_raw_path(
+                year=f"{d.year}", month=f"{d.month:02d}", day=f"{d.day:02d}"
+            )
+            for d in date_list
+        ]
+        filepath_list.sort()
+
+        return filepath_list
+
+    def _get_to_be_loaded_path_list(self):
+        """Get list of filepaths of files to be loaded."""
+        date_list = self._create_date_list()
+
+        filepath_list = [
+            self._get_processed_path(
+                self._get_raw_path(
+                    year=f"{d.year}",
+                    month=f"{d.month:02d}",
+                    day=f"{d.day:02d}",
+                )
+            )
+            for d in date_list
+        ]
+        filepath_list.sort()
+
+        return filepath_list
+
+    @staticmethod
+    def _read_csv_from_url(url):
+
+        context = ssl.create_default_context()
+        context.set_ciphers("DEFAULT")
+        result = urlopen(url, context=context)
+
+        return pd.read_csv(result)
+
+    @check_file_existence
+    def _download(self, filepath: Path, url: str, clobber: bool) -> Path:
+        logger.info("Downloading CHIRPS NetCDF file.")
+        response = requests.get(
+            url,
+        )
+        with open(filepath, "wb") as out_file:
+            out_file.write(response.content)
+        return filepath
+
+    @check_file_existence
+    def _process(self, filepath: Path, ds, clobber: bool) -> Path:
+        pass
+
+
+class ChirpsMonthly(_Chirps):
+    """
+    Class object to retrieve CHIRPS observational monthly precipitation data.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    geo_bounding_box: GeoBoundingBox
+        the bounding coordinates of the area that should be included in the
+        data.
+    start_date: Optional[Union[datetime.date, str]], default = None
+        Data will be considered starting from date `start_date`.
+        Input can be an ISO8601 string or `datetime.date` object.
+        If None, it is set to 1981-1-1.
+    end_date: Optional[Union[datetime.date, str]], default = None
+        Data will be considered up to date `end_date`.
+        Input can be an ISO8601 string or `datetime.date` object.
+        If None, it is set to the date for which most recent data
+        is available.
+
+    Examples
+    --------
+    >>> from aatoolbox import create_country_config, CodAB, GeoBoundingBox
+    >>> from aatoolbox import ChirpsMonthly
+    >>> import datetime
+    >>>
+    >>> country_config = create_country_config(iso3="bfa")
+    >>> codab = CodAB(country_config=country_config)
+    >>> codab.download()
+    >>> admin0 = codab.load(admin_level=0)
+    >>> geo_bounding_box = GeoBoundingBox.from_shape(admin0)
+    >>>
+    >>> start_date = datetime.date(year=2007, month=10, day=23)
+    >>> end_date = datetime.date(year=2020, month=3, day=2)
+    >>> chirps_monthly = ChirpsMonthly(
+    ...   country_config=country_config,
+    ...   geo_bounding_box=geo_bounding_box,
+    ...   start_date=start_date,
+    ...   end_date=end_date
+    ... )
+    >>> chirps_monthly.download()
+    >>> chirps_monthly.process()
+    >>> chirps_monthly_data = chirps_monthly.load()
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        geo_bounding_box: GeoBoundingBox,
+        start_date: Optional[Union[date, str]] = None,
+        end_date: Optional[Union[date, str]] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            geo_bounding_box=geo_bounding_box,
+            date_range_freq="SMS",
+            frequency="monthly",
+            resolution=0.05,
+            start_date=start_date,
+            end_date=end_date,
+        )
+
+    def _get_file_name(
+        self,
+        year: str,
+        month: str,
+        day: str,
+    ) -> str:
+
+        file_name_base = self._get_file_name_base(
+            year=year,
+            month=month,
+        )
+
+        file_name = (
+            f"{file_name_base}"
+            f"r{self._resolution}"
+            f"_{self._geobb.get_filename_repr(p=0)}.nc"
+        )
+        return file_name
+
+    def _get_last_available_date(self):
+        """Get the most recent date for which data is available."""
+        # The url contains a table where the last date available is
+        # specified.
+        url = (
+            f"{_BASE_URL}"
+            ".monthly/.global/"
+            ".T/last/subgrid/0./add/T/table%3A/1/%3Atable/.csv"
+        )
+
+        df = self._read_csv_from_url(url)
+
+        datetime_object = datetime.strptime(df.values[0][0], "%b %Y")
+        day = calendar.monthrange(datetime_object.year, datetime_object.month)[
+            1
+        ]
+        datetime_object = datetime_object.replace(day=day)
+
+        return datetime_object.date()
+
+    def _get_url(self, year: str, month: str, day: str) -> str:
+
+        # Convert month from month number (in string format) to
+        # three-letter name
+        month_name = calendar.month_abbr[int(month)]
+
+        location_url = self._get_location_url()
+
+        url = (
+            f"{_BASE_URL}"
+            ".monthly/.global/.precipitation/"
+            f"{location_url}"
+            f"T/%28{month_name}%20{year}%29%28{month_name}%20{year}"
+            "%29RANGEEDGES/data.nc"
+        )
+
+        return url
+
+    @check_file_existence
+    def _process(self, filepath: Path, ds, clobber: bool) -> Path:
+        # fix dates
+        ds.aat.correct_calendar(inplace=True)
+        ds = xr.decode_cf(ds)
+        if "prcp" in list(ds.keys()):
+            ds = ds.rename({"prcp": "precipitation"})
+        xr.Dataset.to_netcdf(ds, path=filepath)
+        return filepath
+
+
+class ChirpsDaily(_Chirps):
+    """
+    Class object to retrieve CHIRPS observational monthly precipitation data.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    geo_bounding_box: GeoBoundingBox
+        the bounding coordinates of the area that should be included in the
+        data.
+    resolution: float, default = 0.05
+        resolution of data to be downloaded. Can be 0.05 or 0.25.
+    start_date: Optional[Union[datetime.date, str]], default = None
+        Data will be considered starting from date `start_date`.
+        Input can be an ISO8601 string or `datetime.date` object.
+        If None, it is set to 1981-1-1.
+    end_date: Optional[Union[datetime.date, str]], default = None
+        Data will be considered up to date `end_date`.
+        Input can be an ISO8601 string or `datetime.date` object.
+        If None, it is set to the date for which most recent data
+        is available.
+
+    Examples
+    --------
+    >>> from aatoolbox import create_country_config, CodAB, GeoBoundingBox,
+    ... ChirpsDaily
+    >>> import datetime
+    >>>
+    >>> country_config = create_country_config(iso3="bfa")
+    >>> codab = CodAB(country_config=country_config)
+    >>> codab.download()
+    >>> admin0 = codab.load(admin_level=0)
+    >>> geo_bounding_box = GeoBoundingBox.from_shape(admin0)
+    >>> start_date = datetime.date(year=2007, month=10, day=23)
+    >>> end_date = datetime.date(year=2020, month=3, day=2)
+    >>> chirps_daily = ChirpsDaily(
+    ...   country_config=country_config,
+    ...   geo_bounding_box=geo_bounding_box,
+    ...   start_date=start_date,
+    ...   end_date=end_date
+    ... )
+    >>> chirps_daily.download()
+    >>> chirps_daily.process()
+    >>> chirps_daily_data = chirps_daily.load()
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        geo_bounding_box: GeoBoundingBox,
+        resolution: float = 0.05,
+        start_date: Optional[Union[date, str]] = None,
+        end_date: Optional[Union[date, str]] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            geo_bounding_box=geo_bounding_box,
+            date_range_freq="D",
+            frequency="daily",
+            resolution=resolution,
+            start_date=start_date,
+            end_date=end_date,
+        )
+
+    def _get_file_name(
+        self,
+        year: str,
+        month: str,
+        day: str,
+    ) -> str:
+
+        file_name_base = self._get_file_name_base(
+            year=year,
+            month=month,
+        )
+
+        if len(day) == 1:
+            day = f"0{day}"
+
+        file_name = (
+            f"{file_name_base}"
+            f"{day}_"
+            f"r{self._resolution}"
+            f"_{self._geobb.get_filename_repr(p=0)}.nc"
+        )
+        return file_name
+
+    def _get_last_available_date(self):
+        """Get the most recent date for which data is available."""
+        # The url contains a table where the last date available is
+        # specified.
+        url = (
+            f"{_BASE_URL}"
+            ".daily-improved/.global/"
+            f"{str(self._resolution).replace('.', 'p')}/"
+            ".T/last/subgrid/0./add/T/table%3A/1/%3Atable/.csv"
+        )
+
+        df = self._read_csv_from_url(url)
+
+        datetime_object = datetime.strptime(df.values[0][0], "%d %b %Y")
+
+        return datetime_object.date()
+
+    def _get_url(self, year: str, month: str, day: str) -> str:
+
+        # Convert month from month number (in string format) to
+        # three-letter name
+        month_name = calendar.month_abbr[int(month)]
+
+        location_url = self._get_location_url()
+
+        url = (
+            f"{_BASE_URL}"
+            ".daily-improved/.global/."
+            f"{str(self._resolution).replace('.', 'p')}/.prcp/"
+            f"{location_url}"
+            f"T/%28{day}%20{month_name}%20{year}%29%28{day}"
+            f"%20{month_name}%20{year}"
+            "%29RANGEEDGES/data.nc"
+        )
+
+        return url
+
+    @check_file_existence
+    def _process(self, filepath: Path, ds, clobber: bool) -> Path:
+        # fix dates
+        ds = ds.assign_coords(
+            T=cftime.datetime.fromordinal(
+                ds.T.values, calendar="standard", has_year_zero=False
+            )
+        )
+        ds.aat.correct_calendar(inplace=True)
+        if "prcp" in list(ds.keys()):
+            ds = ds.rename({"prcp": "precipitation"})
+        xr.Dataset.to_netcdf(ds, path=filepath)
+        return filepath
diff --git a/src/aatoolbox/datasources/codab/codab.py b/src/aatoolbox/datasources/codab/codab.py
index f90bb459..78b2408a 100644
--- a/src/aatoolbox/datasources/codab/codab.py
+++ b/src/aatoolbox/datasources/codab/codab.py
@@ -108,8 +108,8 @@ def load(self, admin_level: int) -> gpd.GeoDataFrame:  # type: ignore
                 f"config file"
             )
         return self._load_admin_layer(
-            layer_name=self._datasource_config.layer_base_name.format(
-                admin_level=admin_level
+            layer_name=getattr(
+                self._datasource_config, f"admin{admin_level}_name"
             )
         )
 
diff --git a/src/aatoolbox/datasources/datasource.py b/src/aatoolbox/datasources/datasource.py
index d3ecef69..cf651fd2 100644
--- a/src/aatoolbox/datasources/datasource.py
+++ b/src/aatoolbox/datasources/datasource.py
@@ -60,16 +60,23 @@ def __init__(
             is_public=is_public, is_raw=False, is_global=is_global_processed
         )
 
+    @staticmethod
     def _config_attribute_name_validator(
-        self, config_datasource_name: str, country_config: CountryConfig
+        config_datasource_name: str, country_config: CountryConfig
     ):
         try:
-            return getattr(country_config, config_datasource_name)
-        except AttributeError as err:
+            datasource_config = getattr(country_config, config_datasource_name)
+        except AttributeError:
+            datasource_config = None
+        # If the datasource is one of the defaults, it's set to None and
+        # thus an attribute error won't be raised. So also need to check for
+        # the case when it is None.
+        if datasource_config is None:
             raise AttributeError(
-                f"{config_datasource_name} needs to be added to the the "
+                f"{config_datasource_name} needs to be added to the "
                 f"config file. See the documentation for more details."
-            ) from err
+            )
+        return datasource_config
 
     def _get_base_dir(
         self,
diff --git a/src/aatoolbox/datasources/fewsnet/fewsnet.py b/src/aatoolbox/datasources/fewsnet/fewsnet.py
index 412ea315..491cb943 100644
--- a/src/aatoolbox/datasources/fewsnet/fewsnet.py
+++ b/src/aatoolbox/datasources/fewsnet/fewsnet.py
@@ -82,13 +82,6 @@ def __init__(
                 f" {self._country_config.iso3}."
             )
 
-        if self._country_config.fewsnet is None:
-            raise AttributeError(
-                "The country configuration file does not contain "
-                "any FEWS NET region name. Please update the config file and "
-                "try again. See the documentation for the valid region names. "
-            )
-
     # mypy will give error Signature of "download" incompatible with supertype
     # "DataSource" due to `pub_year` and `pub_month` not being an arg in
     # `DataSource`. This is however valid so ignore mypy
diff --git a/src/aatoolbox/datasources/glofas/__init__.py b/src/aatoolbox/datasources/glofas/__init__.py
new file mode 100644
index 00000000..0b0c4207
--- /dev/null
+++ b/src/aatoolbox/datasources/glofas/__init__.py
@@ -0,0 +1,1 @@
+"""Methods pertaining to GloFAS data."""
diff --git a/src/aatoolbox/datasources/glofas/forecast.py b/src/aatoolbox/datasources/glofas/forecast.py
new file mode 100644
index 00000000..a9f2f439
--- /dev/null
+++ b/src/aatoolbox/datasources/glofas/forecast.py
@@ -0,0 +1,269 @@
+"""Download and process GloFAS forecast and reforecast river discharge data."""
+import logging
+from datetime import date
+from pathlib import Path
+from typing import List, Union
+
+import xarray as xr
+from dateutil import rrule
+
+from aatoolbox.config.countryconfig import CountryConfig
+from aatoolbox.datasources.glofas import glofas
+from aatoolbox.utils.check_file_existence import check_file_existence
+from aatoolbox.utils.geoboundingbox import GeoBoundingBox
+
+logger = logging.getLogger(__name__)
+
+
+class _GlofasForecastBase(glofas.Glofas):
+    """Base class for all GloFAS forecast data downloading and processing."""
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        geo_bounding_box: GeoBoundingBox,
+        cds_name: str,
+        system_version: str,
+        product_type: Union[str, List[str]],
+        date_variable_prefix: str,
+        frequency: int,
+        coord_names: List[str],
+        leadtime_max: int,
+        start_date_min: date,
+        end_date_max: date = None,
+        start_date: Union[date, str] = None,
+        end_date: Union[date, str] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            geo_bounding_box=geo_bounding_box,
+            cds_name=cds_name,
+            system_version=system_version,
+            product_type=product_type,
+            date_variable_prefix=date_variable_prefix,
+            frequency=frequency,
+            coord_names=coord_names,
+            start_date_min=start_date_min,
+            end_date_max=end_date_max,
+            start_date=start_date,
+            end_date=end_date,
+            leadtime_max=leadtime_max,
+        )
+
+    @check_file_existence
+    def _load_single_file(
+        self, input_filepath: Path, filepath: Path, clobber: bool
+    ) -> xr.Dataset:
+        return _read_in_ensemble_and_perturbed_datasets(
+            filepath=input_filepath
+        )
+
+
+def _read_in_ensemble_and_perturbed_datasets(filepath: Path):
+    """Read in forecast and reforecast data.
+
+    The GloFAS forecast and reforecast data GRIB files contain two
+    separate datasets: the control member, generated from the most accurate
+    estimate of current conditions, and the perturbed forecast, which
+    contains N ensemble members created by perturbing the control forecast.
+
+    This function reads in both datasets and creates an N+1 (perturbed
+    + control) ensemble.
+    See `this paper <https://hess.copernicus.org/preprints/hess-2020-532/>`__
+    for more details.
+    """
+    ds_list = []
+    for data_type in ["cf", "pf"]:
+        ds = xr.load_dataset(
+            filepath,
+            engine="cfgrib",
+            backend_kwargs={
+                "indexpath": "",
+                "filter_by_keys": {"dataType": data_type},
+            },
+        )
+        # Extra processing require for control forecast
+        if data_type == "cf":
+            ds = ds.expand_dims(dim="number")
+        ds_list.append(ds)
+    return xr.combine_by_coords(ds_list, combine_attrs="drop_conflicts")
+
+
+class GlofasForecast(_GlofasForecastBase):
+    """
+    Class for downloading and processing GloFAS forecast data.
+
+    The GloFAS forecast dataset is a global raster presenting river
+    discharge forecast from 26 May 2021 until present day (updated daily), see
+    `this paper <https://hess.copernicus.org/preprints/hess-2020-532/>`__
+    for more details. While CDS does have version 3 pre-release data
+    from 2020-2021,
+    we understand that there were some small issues that were fixed
+    in the final version, so at this point in time this module
+    does not support downloading the pre-release data.
+
+    This class downloads the raw raster data
+    `from CDS
+    <https://cds.climate.copernicus.eu/cdsapp#!/dataset/cems-glofas-forecast?tab=overview>`__,
+    and processes it from a raster to a datasets of reporting points from the
+    `GloFAS interface
+    <https://www.globalfloods.eu/glofas-forecasting/>`_.
+    Due to the CDS request size limits, separate files are downloaded per
+    day (that contain all requested lead times).
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    geo_bounding_box: GeoBoundingBox
+        The bounding coordinates of the area that should be included
+    leadtime_max: int
+        The maximum desired lead time D in days. All forecast data for lead
+        times 1 to D days are downloaded
+    start_date : Union[date, str], default: date(year=2021, month=5, day=26)
+        The starting date for the dataset. If left blank, defaults to the
+        earliest available date
+    end_date : Union[date, str], default: date.today()
+        The ending date for the dataset. If left blank, defaults to
+        the current date
+    Examples
+    --------
+    Download, process and load GloFAS forecast data for the past month,
+    for a lead time of 15 days.
+
+    >>> from datetime import date
+    >>> from aatoolbox import create_country_config, CodAB, GeoBoundingBox,
+    ... GlofasForecast
+    >>>
+    >>> country_config = create_country_config(iso3="npl")
+    >>> codab = CodAB(country_config=country_config)
+    >>> codab.download()
+    >>> admin_npl = codab.load()
+    >>> geo_bounding_box = GeoBoundingBox.from_shape(admin_npl)
+    >>>
+    >>> glofas_forecast = GlofasForecast(
+    ...     country_config=country_config,
+    ...     geo_bounding_box=geo_bounding_box,
+    ...     leadtime_max=15,
+    ...     end_date=date(year=2022, month=10, day=22),
+    ...     start_date=date(year=2022, month=9, day=22)
+    ... )
+    >>> glofas_forecast.download()
+    >>> glofas_forecast.process()
+    >>>
+    >>> npl_glofas_forecast_reporting_points = glofas_forecast.load()
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        geo_bounding_box: GeoBoundingBox,
+        leadtime_max: int,
+        start_date: Union[date, str] = None,
+        end_date: Union[date, str] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            geo_bounding_box=geo_bounding_box,
+            start_date=start_date,
+            end_date=end_date,
+            cds_name="cems-glofas-forecast",
+            system_version="operational",
+            product_type=["control_forecast", "ensemble_perturbed_forecasts"],
+            date_variable_prefix="",
+            frequency=rrule.DAILY,
+            coord_names=["number", "step"],
+            leadtime_max=leadtime_max,
+            start_date_min=date(year=2021, month=5, day=26),
+        )
+
+    @staticmethod
+    def _preprocess_load(ds: xr.Dataset) -> xr.Dataset:
+        return ds.expand_dims("time")
+
+
+class GlofasReforecast(_GlofasForecastBase):
+    """
+    Class for downloading and processing GloFAS reforecast data.
+
+    The GloFAS reforecast dataset is a global raster presenting river
+    discharnge forecasted from 1999 until 2018, see
+    `this paper <https://hess.copernicus.org/preprints/hess-2020-532/>`_
+    for more details.
+
+    This class downloads the raw raster data
+    `from CDS
+    <https://cds.climate.copernicus.eu/cdsapp#!/dataset/cems-glofas-reforecast?tab=overview>`__,
+    and processes it from a raster to a datasets of reporting points from the
+    `GloFAS interface
+    <https://www.globalfloods.eu/glofas-forecasting/>`_.
+     Due to the CDS request size limits, separate files are downloaded per
+     month (that contain all requested lead times).
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    geo_bounding_box: GeoBoundingBox
+        The bounding coordinates of the area that should be included
+    leadtime_max: int
+        The maximum desired lead time D in days. All forecast data for lead
+        times 1 to D days are downloaded
+    start_date : Union[date, str], default: date(year=1999, month=1, day=1)
+        The starting date for the dataset. If left blank, defaults to the
+        earliest available date
+    end_date : Union[date, str], default: date(year=2018, month=12, day=31)
+        The ending date for the dataset. If left blank, defaults to the
+        last available date
+
+    Examples
+    --------
+    Download, process and load all available GloFAS reforecast data
+    for a lead time of 15 days.
+
+    >>> from aatoolbox import create_country_config, CodAB, GeoBoundingBox,
+    ... GlofasReforecast
+    >>>
+    >>> country_config = create_country_config(iso3="npl")
+    >>> codab = CodAB(country_config=country_config)
+    >>> codab.download()
+    >>> admin_npl = codab.load()
+    >>> geo_bounding_box = GeoBoundingBox.from_shape(admin_npl)
+    >>>
+    >>> glofas_reforecast = GlofasReforecast(
+    ...     country_config=country_config,
+    ...     geo_bounding_box=geo_bounding_box,
+    ...     leadtime_max=15
+    ... )
+    >>> glofas_reforecast.download()
+    >>> glofas_reforecast.process()
+    >>>
+    >>> npl_glofas_reforecast_reporting_points = glofas_reforecast.load()
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        geo_bounding_box: GeoBoundingBox,
+        leadtime_max: int,
+        start_date: Union[date, str] = None,
+        end_date: Union[date, str] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            geo_bounding_box=geo_bounding_box,
+            start_date=start_date,
+            end_date=end_date,
+            cds_name="cems-glofas-reforecast",
+            system_version="version_3_1",
+            product_type=[
+                "control_reforecast",
+                "ensemble_perturbed_reforecasts",
+            ],
+            date_variable_prefix="h",
+            frequency=rrule.MONTHLY,
+            coord_names=["number", "time", "step"],
+            leadtime_max=leadtime_max,
+            start_date_min=date(year=1999, month=1, day=1),
+            end_date_max=date(year=2018, month=12, day=31),
+        )
diff --git a/src/aatoolbox/datasources/glofas/glofas.py b/src/aatoolbox/datasources/glofas/glofas.py
new file mode 100644
index 00000000..93210d13
--- /dev/null
+++ b/src/aatoolbox/datasources/glofas/glofas.py
@@ -0,0 +1,505 @@
+"""Base class for downloading and processing GloFAS river discharge data."""
+import logging
+import time
+from abc import abstractmethod
+from dataclasses import dataclass
+from datetime import date
+from pathlib import Path
+from typing import List, Tuple, Union
+
+import cdsapi
+import numpy as np
+import xarray as xr
+from dateutil import rrule
+
+from aatoolbox.config.countryconfig import CountryConfig
+from aatoolbox.datasources.datasource import DataSource
+from aatoolbox.utils.dates import get_date_from_user_input
+from aatoolbox.utils.geoboundingbox import GeoBoundingBox
+
+_MODULE_BASENAME = "glofas"
+_HYDROLOGICAL_MODEL = "lisflood"
+_RIVER_DISCHARGE_VAR = "dis24"
+_REQUEST_SLEEP_TIME = 60  # seconds
+
+logger = logging.getLogger(__name__)
+
+
+@dataclass
+class _QueryParams:
+    """
+    Class to keep track of CDS query input and output.
+
+    Parameters
+    ----------
+    filepath: Path
+        Full filepath of downloaded CDS file
+    query: dict
+        Output of _get_query() method, to be submitted to the CDS API
+    request_id: str, default = None
+        Request ID from CDS, only set after request has been made
+    downloaded: bool, default = False
+        Whether or not the file has yet been downloaded
+    """
+
+    filepath: Path
+    query: dict
+    request_id: str = None  # type: ignore # TODO: fix
+    downloaded: bool = False
+
+
+class Glofas(DataSource):
+    """
+    Base class for all GloFAS data downloading and processing.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    geo_bounding_box: GeoBoundingBox
+        The bounding coordinates of the area that should be included
+    cds_name : str
+        The name of the dataset in CDS
+    system_version : str
+        An input to the CDS query
+    product_type : str or list
+        Which product types from the dataset are requested
+    date_variable_prefix : str
+        Some dates require a prefix for the CDS API query
+    frequency : str
+        How to split the query (and thus files): in years, months, or days.
+        Depends on the maximum query size of the product
+    coord_names : list
+        Coordinate names in the xarray dataset
+    start_date_min: date
+        The minimum allowed start date
+    end_date_max: date, default = None
+        The maximum allowed end date
+    start_date : Union[date, str], default = None
+        The starting date for the dataset
+    end_date : Union[date, str], default = None
+        The ending date for the dataset
+    leadtime_max : int, default = None
+        The maximum lead time in days, for forecast or reforecast data
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        geo_bounding_box: GeoBoundingBox,
+        cds_name: str,
+        system_version: str,
+        product_type: Union[str, List[str]],
+        date_variable_prefix: str,
+        frequency: int,
+        coord_names: List[str],
+        start_date_min: date,
+        end_date_max: date = None,
+        start_date: Union[date, str] = None,
+        end_date: Union[date, str] = None,
+        leadtime_max: int = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            datasource_base_dir=_MODULE_BASENAME,
+            is_public=True,
+        )
+        # The GloFAS API on CDS requires coordinates have the format x.x5
+        self._geo_bounding_box = geo_bounding_box.round_coords(
+            offset_val=0.05, round_val=0.1
+        )
+        self._start_date, self._end_date = _set_dates(
+            start_date_min=start_date_min,
+            end_date_max=end_date_max,
+            start_date=start_date,
+            end_date=end_date,
+        )
+        self._cds_name = cds_name
+        self._system_version = system_version
+        self._product_type = product_type
+        self._date_variable_prefix = date_variable_prefix
+        self._frequency = frequency
+        self._coord_names = coord_names
+        self._leadtime_max = leadtime_max
+        self._forecast_type = type(self).__name__.replace("Glofas", "").lower()
+        self._date_range = rrule.rrule(
+            freq=self._frequency,
+            dtstart=self._start_date,
+            until=self._end_date,
+        )
+
+    def download(  # type: ignore
+        self,
+        clobber: bool = False,
+    ) -> List[Path]:
+        """
+        Download the GloFAS data by querying CDS.
+
+        The raw GloFAS data is available as a global raster in CDS. This method
+        downloads the raster files for the specified region of interest
+        and date range. The files are in GRIB format and are split up either
+        by day, month, or year depending on the GloFAS product.
+
+        Parameters
+        ----------
+        clobber : bool, default = False
+            Overwrite files that were already downloaded
+
+        Returns
+        -------
+        A list paths of downloaded files
+
+        """
+        msg = (
+            f"Downloading GloFAS {self._forecast_type} "
+            f"for {self._start_date} - {self._end_date}"
+        )
+        if self._leadtime_max is not None:
+            msg += f"and up to {self._leadtime_max} day lead time"
+        logger.info(msg)
+
+        # Make directory
+        output_directory = self._get_directory()
+        output_directory.mkdir(parents=True, exist_ok=True)
+
+        # Get list of files to open
+        query_params_list = []
+        for file_date in self._date_range:
+            output_filepath = self._get_filepath(
+                year=file_date.year,
+                month=file_date.month,
+                day=file_date.day,
+            )
+            if not clobber and output_filepath.exists():
+                continue
+            query_params_list.append(
+                _QueryParams(
+                    filepath=output_filepath,
+                    query=self._get_query(
+                        year=file_date.year,
+                        month=file_date.month,
+                        day=file_date.day,
+                    ),
+                )
+            )
+
+        download_filepaths = self._download(
+            query_params_list=query_params_list
+        )
+        logger.info(
+            f"Downloaded {len(download_filepaths)} files to {output_directory}"
+        )
+        logger.debug(f"Files downloaded: {download_filepaths}")
+        return download_filepaths
+
+    def process(  # type: ignore
+        self,
+        clobber: bool = False,
+    ) -> List[Path]:
+        """
+        Process the downloaded GloFAS files.
+
+        For each raw GRIB file, read it in and extract the river discharge
+        from the reporting point coordinates specified in the configuration
+        file. Saves the output as a NetCDF file, where files are split by day,
+        month or year depending on the GloFAS product.
+
+        Parameters
+        ----------
+        clobber : bool, default = False
+            Overwrite files that were already processed
+
+        Returns
+        -------
+        A list paths of processed files
+
+        """
+        logger.info(
+            f"Processing GloFAS {self._forecast_type} for "
+            f"{self._start_date} - {self._end_date} and up to "
+            f"{self._leadtime_max} day lead time"
+        )
+        # Make the directory
+        output_directory = self._get_directory(is_processed=True)
+        output_directory.mkdir(parents=True, exist_ok=True)
+        # Get list of files to open
+        processed_filepaths = []
+        for file_date in self._date_range:
+            input_filepath = self._get_filepath(
+                year=file_date.year,
+                month=file_date.month,
+                day=file_date.day,
+            )
+            output_filepath = self._get_filepath(
+                year=file_date.year,
+                month=file_date.month,
+                day=file_date.day,
+                is_processed=True,
+            )
+            if not clobber and output_filepath.exists():
+                continue
+            logger.debug(f"Processing {input_filepath}")
+            ds_raw = self._load_single_file(
+                input_filepath=input_filepath,
+                filepath=output_filepath,
+                clobber=clobber,
+            )
+            ds_processed = self._get_reporting_point_dataset(ds=ds_raw)
+            # NetCDF doesn't like to overwrite files
+            if output_filepath.exists():
+                output_filepath.unlink()
+            ds_processed.to_netcdf(output_filepath)
+            processed_filepaths.append(output_filepath)
+            logger.debug(f"Wrote file to {output_filepath}")
+        logger.info(
+            f"Processed {len(processed_filepaths)} files to {output_directory}"
+        )
+        logger.debug(f"Files downloaded: {processed_filepaths}")
+        return processed_filepaths
+
+    def load(
+        self,
+    ) -> xr.Dataset:
+        """
+        Load the processed GloFAS data as an xarray.DataSet.
+
+        Returns
+        -------
+        A single xarray dataset containing all GloFAS reporting points
+        and their associated river discharge
+
+        """
+        filepath_list = [
+            self._get_filepath(
+                year=dataset_date.year,
+                month=dataset_date.month,
+                day=dataset_date.day,
+                is_processed=True,
+            )
+            for dataset_date in self._date_range
+        ]
+        with xr.open_mfdataset(
+            filepath_list, preprocess=self._preprocess_load
+        ) as ds:
+            return ds
+
+    def _get_filepath(
+        self,
+        year: int,
+        month: int = None,
+        day: int = None,
+        is_processed: bool = False,
+    ) -> Path:
+        """Get downloaded / processed filepaths based on GloFAS product."""
+        filename = f"{self._country_config.iso3}_{self._cds_name}_{year}"
+        if self._frequency in [rrule.MONTHLY, rrule.DAILY]:
+            filename += f"-{str(month).zfill(2)}"
+        if self._frequency == rrule.DAILY:
+            filename += f"-{str(day).zfill(2)}"
+        if self._leadtime_max is not None:
+            filename += f"_ltmax{str(self._leadtime_max).zfill(2)}d"
+        filename += f"_{self._geo_bounding_box.get_filename_repr(p=2)}"
+        if is_processed:
+            filename += "_processed.nc"
+        else:
+            filename += ".grib"
+        return self._get_directory(is_processed=is_processed) / Path(filename)
+
+    def _get_directory(self, is_processed: bool = False) -> Path:
+        """Get download / processed directory for GloFAS product."""
+        return (
+            self._processed_base_dir
+            if is_processed
+            else self._raw_base_dir / self._cds_name
+        )
+
+    def _download(
+        self,
+        query_params_list: List[_QueryParams],
+    ) -> List[Path]:
+        """
+        Download the GloFAS data from CDS.
+
+        Uses query_params_list, which is a list of API request input dicts,
+        to query the CDS API, and for each query,  store the request
+        ID that is returned, and the downloaded state.
+
+        Then loops through the list of request, checking each one to see
+        if it has been completed on the CDS side. If so, it's downloaded,
+        and then removed from the list of requests. The process continues
+        until the request list is empty
+        """
+        # First make the requests to the CDS client and store request number
+        for query_params in query_params_list:
+            logger.debug(f"Making request {query_params.query}")
+            query_params.request_id = (
+                cdsapi.Client(wait_until_complete=False, delete=False)
+                .retrieve(name=self._cds_name, request=query_params.query)
+                .reply["request_id"]
+            )
+        # Loop through the request list and check status until all requests
+        # are downloaded
+        downloaded_filepaths = []
+        while query_params_list:
+            for query_params in query_params_list:
+                result = cdsapi.api.Result(
+                    client=cdsapi.Client(
+                        wait_until_complete=False, delete=False
+                    ),
+                    reply={"request_id": query_params.request_id},
+                )
+                result.update()
+                state = result.reply["state"]
+                logger.debug(
+                    f"For request {query_params.request_id} and filename "
+                    f"{query_params.filepath}, state is {state}"
+                )
+                if state == "completed":
+                    result.download(query_params.filepath)
+                    query_params.downloaded = True
+                    downloaded_filepaths.append(query_params.filepath)
+                elif state == "failed":
+                    raise RuntimeError("Query has failed, try again")
+            # Remove requests that have been downloaded
+            query_params_list = [
+                query_params
+                for query_params in query_params_list
+                if not query_params.downloaded
+            ]
+            # Sleep a bit before the next loop so that we're not
+            # hammering on cds
+            if query_params_list:
+                time.sleep(_REQUEST_SLEEP_TIME)
+        return downloaded_filepaths
+
+    def _get_query(
+        self,
+        year: int,
+        month: int = None,
+        day: int = None,
+    ) -> dict:
+        """Create dictionary for CDS API query input."""
+        query = {
+            "variable": "river_discharge_in_the_last_24_hours",
+            "format": "grib",
+            "product_type": self._product_type,
+            "system_version": self._system_version,
+            "hydrological_model": _HYDROLOGICAL_MODEL,
+            f"{self._date_variable_prefix}year": str(year),
+            f"{self._date_variable_prefix}month": str(month).zfill(2)
+            if self._frequency in [rrule.MONTHLY, rrule.DAILY]
+            else [str(x + 1).zfill(2) for x in range(12)],
+            f"{self._date_variable_prefix}day": str(day).zfill(2)
+            if self._frequency == rrule.DAILY
+            else [str(x + 1).zfill(2) for x in range(31)],
+            "area": [
+                self._geo_bounding_box.lat_max,
+                self._geo_bounding_box.lon_min,
+                self._geo_bounding_box.lat_min,
+                self._geo_bounding_box.lon_max,
+            ],
+        }
+        if self._leadtime_max is not None:
+            leadtime = list(np.arange(self._leadtime_max) + 1)
+            query["leadtime_hour"] = [
+                str(single_leadtime * 24) for single_leadtime in leadtime
+            ]
+        logger.debug(f"Query: {query}")
+        return query
+
+    @abstractmethod
+    def _load_single_file(self, *args, **kwargs) -> xr.Dataset:
+        """Process a single raw raster file."""
+        pass
+
+    def _get_reporting_point_dataset(self, ds: xr.Dataset) -> xr.Dataset:
+        """Convert raw raster to processed that uses reporting points."""
+        if self._country_config.glofas is None:
+            raise KeyError(
+                "The country configuration file does not contain "
+                "any reporting point coordinates. Please update the "
+                "configuration file and try again."
+            )
+        # Check that lat and lon of reporting points are in the bounds
+        for reporting_point in self._country_config.glofas.reporting_points:
+            if (
+                not ds.longitude.min()
+                < reporting_point.lon
+                < ds.longitude.max()
+            ):
+                raise IndexError(
+                    f"ReportingPoint {reporting_point.id} has out-of-bounds "
+                    f"lon value of {reporting_point.lon} (data lon ranges "
+                    f"from {ds.longitude.min().values} "
+                    f"to {ds.longitude.max().values})"
+                )
+            if not ds.latitude.min() < reporting_point.lat < ds.latitude.max():
+                raise IndexError(
+                    f"ReportingPoint {reporting_point.id} has out-of-bounds "
+                    f"lat value of {reporting_point.lat} (data lat ranges "
+                    f"from {ds.latitude.min().values} "
+                    f"to {ds.latitude.max().values})"
+                )
+        # If reporting points fit then return processed dataset
+        return xr.Dataset(
+            data_vars={
+                reporting_point.name: (
+                    self._coord_names,
+                    ds.sel(
+                        longitude=reporting_point.lon,
+                        latitude=reporting_point.lat,
+                        method="nearest",
+                    )[_RIVER_DISCHARGE_VAR].data,
+                )
+                # fmt: off
+                for reporting_point in
+                self._country_config.glofas.reporting_points
+                # fmt: on
+            },
+            coords={
+                coord_name: ds[coord_name] for coord_name in self._coord_names
+            },
+        )
+
+    @staticmethod
+    def _preprocess_load(ds: xr.Dataset) -> xr.Dataset:
+        """Preprocessing to do before loading the processed data."""
+        return ds
+
+
+def _set_dates(
+    start_date_min: date,
+    end_date_max: date = None,
+    start_date: Union[date, str] = None,
+    end_date: Union[date, str] = None,
+) -> Tuple[date, date]:
+    """Adjust date types and check against limits."""
+    # TODO: perhaps this is general enough to be useful for other data sources?
+    # Set any missing dates
+    if end_date_max is None:
+        end_date_max = date.today()
+    if start_date is None:
+        start_date = start_date_min
+    if end_date is None:
+        end_date = end_date_max
+    # Make sure that dates are all date type
+    start_date = get_date_from_user_input(start_date)
+    end_date = get_date_from_user_input(end_date)
+    # Check against bounds
+    if start_date > end_date:
+        raise ValueError(
+            "Please ensure that the start date is <= the end_date"
+        )
+    if start_date < start_date_min:
+        logger.warning(
+            f"Start date {start_date} is too far in the past,"
+            f"setting to {start_date_min}"
+        )
+        start_date = start_date_min
+    if end_date > end_date_max:
+        logger.warning(
+            f"End date {end_date} is too far in the future,"
+            f"setting to {end_date_max}"
+        )
+        end_date = end_date_max
+    return start_date, end_date
diff --git a/src/aatoolbox/datasources/glofas/reanalysis.py b/src/aatoolbox/datasources/glofas/reanalysis.py
new file mode 100644
index 00000000..a4636d4a
--- /dev/null
+++ b/src/aatoolbox/datasources/glofas/reanalysis.py
@@ -0,0 +1,101 @@
+"""Download and process GloFAS reanalysis river discharge data."""
+import logging
+from datetime import date
+from pathlib import Path
+from typing import Union
+
+import xarray as xr
+from dateutil import rrule
+
+from aatoolbox.config.countryconfig import CountryConfig
+from aatoolbox.datasources.glofas import glofas
+from aatoolbox.utils.check_file_existence import check_file_existence
+from aatoolbox.utils.geoboundingbox import GeoBoundingBox
+
+logger = logging.getLogger(__name__)
+
+
+class GlofasReanalysis(glofas.Glofas):
+    """
+    Class for downloading and processing GloFAS reanalysis data.
+
+    The GloFAS reanalysis dataset is a global raster presenting river
+    discharnge from 1979 until present day (updated daily), see
+    `this paper <https://essd.copernicus.org/articles/12/2043/2020/>`_
+    for more details.
+
+    This class downloads the raw raster data
+    `from CDS
+    <https://cds.climate.copernicus.eu/cdsapp#!/dataset/cems-glofas-historical?tab=overview>`_,
+    and processes it from a raster to a datasets of reporting points from the
+    `GloFAS interface
+    <https://www.globalfloods.eu/glofas-forecasting/>`_.
+    Due to the CDS request size limits, separate files are downloaded per
+    year.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    geo_bounding_box: GeoBoundingBox
+        The bounding coordinates of the area that should be included
+    start_date : Union[date, str], default: date(1979, 1, 1)
+        The starting date for the dataset. If left blank, defaults to the
+        earliest available date
+    end_date : Union[date, str], default: date.today()
+        The ending date for the dataset. If left blank, defaults to
+        the current date
+    Examples
+    --------
+    Download, process and load all historical GloFAS reanalysis data
+    until the current date, set to Oct 22, 2022 for this example.
+
+    >>> from datetime import date
+    >>> from aatoolbox import create_country_config, CodAB, GeoBoundingBox,
+    ... GlofasReanalysis
+    >>>
+    >>> country_config = create_country_config(iso3="bgd")
+    >>> codab = CodAB(country_config=country_config)
+    >>> codab.download()
+    >>> admin_npl = codab.load()
+    >>> geo_bounding_box = GeoBoundingBox.from_shape(admin_npl)
+    >>>
+    >>> glofas_reanalysis = GlofasReanalysis(
+    ...     country_config=country_config,
+    ...     geo_bounding_box=geo_bounding_box,
+    ...     end_date=date(year=2022, month=10, day=22)
+    ... )
+    >>> glofas_reanalysis.download()
+    >>> glofas_reanalysis.process()
+    >>>
+    >>> npl_glofas_reanalysis_reporting_points = glofas_reanalysis.load()
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        geo_bounding_box: GeoBoundingBox,
+        start_date: Union[date, str] = None,
+        end_date: Union[date, str] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            geo_bounding_box=geo_bounding_box,
+            cds_name="cems-glofas-historical",
+            system_version="version_3_1",
+            product_type="consolidated",
+            date_variable_prefix="h",
+            frequency=rrule.YEARLY,
+            coord_names=["time"],
+            start_date_min=date(year=1979, month=1, day=1),
+            start_date=start_date,
+            end_date=end_date,
+        )
+
+    @check_file_existence
+    def _load_single_file(
+        self, input_filepath: Path, filepath: Path, clobber: bool
+    ) -> xr.Dataset:
+        return xr.load_dataset(
+            input_filepath, engine="cfgrib", backend_kwargs={"indexpath": ""}
+        )
diff --git a/src/aatoolbox/datasources/iri/iri_seasonal_forecast.py b/src/aatoolbox/datasources/iri/iri_seasonal_forecast.py
index 232c6d7e..48893e9c 100644
--- a/src/aatoolbox/datasources/iri/iri_seasonal_forecast.py
+++ b/src/aatoolbox/datasources/iri/iri_seasonal_forecast.py
@@ -65,7 +65,7 @@ def download(
         To download data from the IRI API, a key is required for
         authentication, and must be set in the ``IRI_AUTH`` environment
         variable. To obtain this key config you need to create an account
-        `here.<https://iridl.ldeo.columbia.edu/auth/login>`_.
+        `here <https://iridl.ldeo.columbia.edu/auth/login>`_.
         Note that this key might be changed over time, and need to be updated
         regularly.
 
@@ -179,15 +179,6 @@ def _load_raw(self) -> xr.Dataset:
                 drop_variables="C",
             )
             return ds
-        except ValueError as err:
-            raise ValueError(
-                f"Cannot open the netcdf file {self._get_raw_path()}. "
-                f"This might be due to invalid "
-                f"download with wrong authentication. Check the validity of "
-                f"the authentication key found in your {_IRI_AUTH} environment"
-                f"variable and try to download again. Otherwise make sure the "
-                f"correct backend for opening a netCDF file is installed."
-            ) from err
         except FileNotFoundError as err:
             raise FileNotFoundError(
                 f"Cannot open the netcdf file {self._get_raw_path()}. Make "
@@ -205,6 +196,15 @@ def _download(
             # have to authenticate by using a cookie
             cookies={"__dlauth_id": iri_auth},
         )
+        if response.headers["Content-Type"] != "application/x-netcdf":
+            msg = (
+                f"The request returned headers indicating that the expected "
+                f"file type was not returned. In some cases th  is may be due "
+                f"to an issue with the authentication. Please check the "
+                f"validity of the authentication key found in your "
+                f"{_IRI_AUTH} environment variable and try again."
+            )
+            raise requests.RequestException(msg)
         with open(filepath, "wb") as out_file:
             out_file.write(response.content)
         return filepath
diff --git a/src/aatoolbox/datasources/usgs/__init__.py b/src/aatoolbox/datasources/usgs/__init__.py
new file mode 100644
index 00000000..4ba5f4cc
--- /dev/null
+++ b/src/aatoolbox/datasources/usgs/__init__.py
@@ -0,0 +1,1 @@
+"""Methods pertaining to NDVI data produced by USGS."""
diff --git a/src/aatoolbox/datasources/usgs/ndvi_base.py b/src/aatoolbox/datasources/usgs/ndvi_base.py
new file mode 100644
index 00000000..059bb201
--- /dev/null
+++ b/src/aatoolbox/datasources/usgs/ndvi_base.py
@@ -0,0 +1,848 @@
+"""Class to download and process USGS eMODIS NDVI data.
+
+Download, process, and load eMODIS NDVI data published
+in the `USGS FEWS NET data portal
+<https://earlywarning.usgs.gov/fews>`_.
+
+Warning: the MODIS sensor has been reported by USGS to
+have degraded in quality since May 2022 (dekad 13), and
+updates to this data source have stopped. This module
+remains for users to have access to historic data but
+recent data is unavailable and care should be used
+analyzing any data since dekad 13 of 2022.
+"""
+
+# TODO: add progress bar
+import functools
+import logging
+import re
+from datetime import date, datetime
+from io import BytesIO
+from pathlib import Path
+from typing import Generator, List, Optional, Tuple, Union
+from urllib.error import HTTPError
+from urllib.request import urlopen
+from zipfile import ZipFile
+
+import geopandas as gpd
+import pandas as pd
+import rioxarray  # noqa: F401
+import xarray as xr
+from rasterio.errors import RasterioIOError
+
+import aatoolbox.utils.raster  # noqa: F401
+from aatoolbox.config.countryconfig import CountryConfig
+from aatoolbox.datasources.datasource import DataSource
+from aatoolbox.utils.dates import (
+    compare_dekads_gt,
+    compare_dekads_lt,
+    dekad_to_date,
+    expand_dekads,
+    get_dekadal_date,
+)
+
+logger = logging.getLogger(__name__)
+
+_DATE_TYPE = Union[date, str, Tuple[int, int], None]
+_EARLIEST_DATE = (2002, 19)
+
+# USGS has reported degradation of USGS NDVI data
+# from the below date and warnings should be used
+_DEGRADATION_DATE = (2022, 13)
+
+
+class _UsgsNdvi(DataSource):
+    """Base class to retrieve USGS NDVI data.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    data_variable : str
+        Data variable date
+    data_variable_suffix : str
+        Data variable file string
+    data_variable_url : str
+        URL string for data variable
+    start_date : _DATE_TYPE, default = None
+        Start date. Can be passed as a ``datetime.date``
+        object or a data string in ISO8601 format, and
+        the relevant dekad will be determined. Or pass
+        directly as year-dekad tuple, e.g. (2020, 1).
+        If ``None``, ``start_date`` is set to earliest
+        date with data: 2002, dekad 19.
+    end_date : _DATE_TYPE, default = None
+        End date. Can be passed as a ``datetime.date``
+        object and the relevant dekad will be determined,
+        as a date string in ISO8601 format, or as a
+        year-dekad tuple, i.e. (2020, 1). If ``None``,
+        ``end_date`` is set to ``date.today()``.
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        data_variable: str,
+        data_variable_suffix: str,
+        data_variable_url: str,
+        start_date: _DATE_TYPE = None,
+        end_date: _DATE_TYPE = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            datasource_base_dir="usgs_ndvi",
+            is_public=True,
+            is_global_raw=True,
+            config_datasource_name="usgs_ndvi",
+        )
+
+        # set data variable
+        self._data_variable = data_variable
+        self._data_variable_url = data_variable_url
+        self._data_variable_suffix = data_variable_suffix
+
+        # set dates for data download and processing
+        self._start_date = get_dekadal_date(
+            input_date=start_date, default_date=_EARLIEST_DATE
+        )
+
+        self._end_date = get_dekadal_date(
+            input_date=end_date, default_date=date.today()
+        )
+
+        if compare_dekads_gt(self._end_date, _DEGRADATION_DATE):
+            logger.warning(
+                "USGS has reported degradation of eMODIS NDVI data "
+                "due to issues with the MODIS sensor's satellite. "
+                "This affects NDVI data from May 2022 (dekad 13), "
+                "and users should be very careful of using any "
+                "results after this date. This module is maintained "
+                "for users to have access to historic data."
+            )
+
+        # warn if dates outside earliest dates
+        if compare_dekads_lt(self._start_date, _EARLIEST_DATE):
+            logger.warning(
+                "Start date is before earliest date data is available. "
+                f"Data will be downloaded from {_EARLIEST_DATE[0]}, dekad "
+                f"{_EARLIEST_DATE[1]}."
+            )
+
+    def download(self, clobber: bool = False) -> Path:
+        """Download raw NDVI data as .tif files.
+
+        NDVI data is downloaded from the USGS API,
+        with data for individual regions, years, and
+        dekads stored as separate .tif files. No
+        authentication is required. Data is downloaded
+        for all available dekads from ``self._start_date``
+        to ``self._end_date``.
+
+        Parameters
+        ----------
+        clobber : bool, default = False
+            If True, overwrites existing files
+
+        Returns
+        -------
+        Path
+            The downloaded filepath
+
+        Examples
+        --------
+        >>> from aatoolbox import create_country_config, \
+        ...  CodAB, UsgsNdviSmoothed
+        >>>
+        >>> # Retrieve admin 2 boundaries for Burkina Faso
+        >>> country_config = create_country_config(iso3="bfa")
+        >>>
+        >>> # setup NDVI
+        >>> bfa_ndvi = UsgsNdviSmoothed(
+        ...     country_config=country_config,
+        ...     start_date=[2020, 1],
+        ...     end_date=[2020, 3]
+        ... )
+        >>> bfa_ndvi.download()
+        """
+        download_dekads = expand_dekads(
+            dekad1=self._start_date, dekad2=self._end_date
+        )
+        for year, dekad in download_dekads:
+            self._download_ndvi_dekad(year=year, dekad=dekad, clobber=clobber)
+        return self._raw_base_dir
+
+    def process(  # type: ignore
+        self,
+        gdf: gpd.GeoDataFrame,
+        feature_col: str,
+        clobber: bool = False,
+        **kwargs,
+    ) -> Path:
+        """
+        Process NDVI data for specific area.
+
+        NDVI data is clipped to the provided
+        ``geometries``, usually a geopandas
+        dataframes ``geometry`` feature. ``kwargs``
+        are passed on to ``aat.computer_raster_stats()``.
+        The ``feature_col`` is used to define
+        the unique processed file.
+
+        The processing keeps track of the latest timestamp of
+        when the raw raster files were modified. If the latest
+        timestamp of the raw data is greater than when it was
+        last processed, the file will automatically be
+        re-processed.
+
+        Parameters
+        ----------
+        gdf : geopandas.GeoDataFrame
+            GeoDataFrame with row per area for stats computation.
+            If ``pd.DataFrame`` is passed, geometry column must
+            have the name ``geometry``. Passed to
+            ``aat.compute_raster_stats()``.
+        feature_col : str
+            Column in ``gdf`` to use as row/feature identifier.
+            and dates. Passed to ``aat.compute_raster_stats()``.
+            The string is also used as a suffix to the
+            processed file path for unique identication of
+            analyses done on different files and columns.
+        clobber : bool, default = False
+            If True, overwrites existing processed dates. If
+            the new file matches the old file, dates will be
+            reprocessed and appended to the data frame. If
+            files do not match, the old file will be replaced.
+            If False, stats are only calculated for year-dekads
+            that have not already been calculated within the
+            file. However, if False and files do not match,
+            value error will be raised.
+        **kwargs
+            Additional keyword arguments passed to
+            ``aat.computer_raster_stats()``.
+
+        Returns
+        -------
+        Path
+            The processed path
+
+        Examples
+        --------
+        >>> from aatoolbox import create_country_config, \
+        ...  CodAB, UsgsNdviSmoothed
+        >>>
+        >>> # Retrieve admin 2 boundaries for Burkina Faso
+        >>> country_config = create_country_config(iso3="bfa")
+        >>> codab = CodAB(country_config=country_config)
+        >>> bfa_admin2 = codab.load(admin_level=2)
+        >>> bfa_admin1 = codab.load(admin_level=1)
+        >>>
+        >>> # setup NDVI
+        >>> bfa_ndvi = UsgsNdviSmoothed(
+        ...     country_config=country_config,
+        ...     start_date=[2020, 1],
+        ...     end_date=[2020, 3]
+        ... )
+        >>> bfa_ndvi.download()
+        >>> bfa_ndvi.process(
+        ...     gdf=bfa_admin2,
+        ...     feature_col="ADM2_FR"
+        ... )
+        >>>
+        >>> # process for admin1
+        >>> bfa_ndvi.process(
+        ...     gdf=bfa_admin1,
+        ...     feature_col="ADM1_FR"
+        ... )
+        """
+        # get statistics for processing
+        stats_list = kwargs.pop(
+            "stats_list", ["mean", "std", "min", "max", "sum", "count"]
+        )
+        percentile_list = kwargs.pop("percentile_list", None)
+
+        if not stats_list and not percentile_list:
+            raise ValueError(
+                "At least one of `stats_list` or "
+                "`percentile_list` must not be `None` "
+                "when passed to `process()`."
+            )
+
+        # getting stats and percentiles for processing
+        stats_list = [] if not stats_list else stats_list
+        percentile_list = (
+            []
+            if not percentile_list
+            else [str(percent) for percent in percentile_list]
+        )
+        process_stats = stats_list + percentile_list
+        percentile_identifier = [False] * len(stats_list) + [True] * len(
+            percentile_list
+        )
+
+        # get dates for processing
+        all_dates_to_process = expand_dekads(
+            dekad1=self._start_date, dekad2=self._end_date
+        )
+
+        for stat, is_percentile in zip(process_stats, percentile_identifier):
+            self._process(
+                clobber=clobber,
+                gdf=gdf,
+                feature_col=feature_col,
+                dates_to_process=all_dates_to_process,
+                stat=stat,
+                is_percentile=is_percentile,
+                kwargs=kwargs,
+            )
+
+        return self._processed_base_dir
+
+    def load(self, feature_col: str) -> pd.DataFrame:  # type: ignore
+        """
+        Load the processed USGS NDVI data.
+
+        Parameters
+        ----------
+        feature_col : str
+            String is  used as a suffix to the
+            processed file path for unique identication of
+            analyses done on different files and columns.
+            The same value must be passed to ``process()``.
+
+        Returns
+        -------
+        pd.DataFrame
+            The processed NDVI dataset.
+
+        Raises
+        ------
+        FileNotFoundError
+            If the requested file cannot be found.
+
+        Examples
+        --------
+        >>> from aatoolbox import create_country_config, \
+        ...  CodAB, UsgsNdviSmoothed
+        >>>
+        >>> # Retrieve admin 2 boundaries for Burkina Faso
+        >>> country_config = create_country_config(iso3="bfa")
+        >>> codab = CodAB(country_config=country_config)
+        >>> bfa_admin2 = codab.load(admin_level=2)
+        >>>
+        >>> # setup NDVI
+        >>> bfa_ndvi = UsgsNdviSmoothed(
+        ...     country_config=country_config,
+        ...     start_date=[2020, 1],
+        ...     end_date=[2020, 3]
+        ... )
+        >>> bfa_ndvi.download()
+        >>> bfa_ndvi.process(
+        ...    gdf=bfa_admin2,
+        ...    feature_col="ADM2_FR"
+        )
+        >>> bfa_ndvi.load(feature_col="ADM2_FR")
+        """
+        try:
+            processed_files = self._find_processed_files(
+                feature_col=feature_col
+            )
+            processed_dfs = [
+                self._load(filepath=fp, drop_modified=True)
+                for fp in processed_files
+            ]
+
+            df = functools.reduce(
+                lambda df1, df2: pd.merge(
+                    df1, df2, on=["date", "year", "dekad", feature_col]
+                ),
+                processed_dfs,
+            )
+        except TypeError as err:
+            raise FileNotFoundError(
+                "Files not found to load. Ensure the download() and process() "
+                "methods are called prior to load()."
+            ) from err
+
+        # filter loaded data frame between our instances dates
+        load_dates = expand_dekads(
+            dekad1=self._start_date, dekad2=self._end_date
+        )
+
+        loaded_dates = df[["year", "dekad"]].values.tolist()
+        keep_rows = [tuple(d) in load_dates for d in loaded_dates]
+
+        df = df.loc[keep_rows]
+
+        # put feature column as 1st column
+        df_feature_col = df.pop(feature_col)
+        df.insert(loc=0, column=feature_col, value=df_feature_col)
+
+        return df
+
+    def load_raster(
+        self, load_date: Union[date, str, Tuple[int, int]]
+    ) -> xr.DataArray:
+        """Load raster for specific year and dekad.
+
+        Parameters
+        ----------
+        load_date : Union[date, str, Tuple[int, int]]
+            Date. Can be passed as a ``datetime.date``
+            object and the relevant dekad will be determined,
+            as a date string in ISO8601 format, or as a
+            year-dekad tuple, i.e. (2020, 1).
+
+        Returns
+        -------
+        xr.DataArray
+            Data array of NDVI data.
+
+        Raises
+        ------
+        FileNotFoundError
+            If the requested file cannot be found.
+        """
+        year, dekad = get_dekadal_date(input_date=load_date)
+
+        filepath = self._get_raw_path(year=year, dekad=dekad, local=True)
+        try:
+            da = rioxarray.open_rasterio(filepath)
+            # get time file was updated
+            file_time = datetime.fromtimestamp(filepath.stat().st_mtime)
+            # assign coordinates for year/dekad
+            # time dimension
+            da = (
+                da.assign_coords(
+                    {
+                        "year": year,
+                        "dekad": dekad,
+                        "date": dekad_to_date(dekad=(year, dekad)),
+                        "modified": file_time,
+                    }
+                )
+                .expand_dims("date")
+                .squeeze("band", drop=True)
+            )
+
+            return da
+
+        except RasterioIOError as err:
+            # check if the requested date is outside the instance bounds
+            # don't prevent loading, but use for meaningful error
+            gt_end = compare_dekads_gt(
+                dekad1=(year, dekad), dekad2=self._end_date
+            )
+            lt_start = compare_dekads_lt(
+                dekad1=(year, dekad), dekad2=self._start_date
+            )
+            if gt_end or lt_start:
+                file_warning = (
+                    f"The requested year and dekad, {year}-{dekad}"
+                    f"are {'greater' if gt_end else 'less'} than the "
+                    f"instance {'end' if gt_end else 'start'} year and dekad, "
+                    f"{self._end_date[0] if gt_end else self._start_date[0]}-"
+                    f"{self._end_date[1] if gt_end else self._start_date[1]}. "
+                    "Calling the `download()` method will not download this "
+                    "file, and you need to re-instantiate the class to "
+                    "include these dates."
+                )
+            else:
+                file_warning = (
+                    "Make sure that you have called the `download()` "
+                    f"method and that the file {filepath.name} exists "
+                    f"in {filepath.parent}."
+                )
+            raise FileNotFoundError(
+                f"Cannot open the .tif file {filepath}. {file_warning}"
+            ) from err
+
+    def _download_ndvi_dekad(
+        self, year: int, dekad: int, clobber: bool
+    ) -> None:
+        """Download NDVI for specific dekad.
+
+        Parameters
+        ----------
+        year : int
+            Year
+        dekad : int
+            Dekad
+        clobber : bool
+            If True, overwrites existing file
+        """
+        filepath = self._get_raw_path(year=year, dekad=dekad, local=True)
+        url_filename = self._get_raw_filename(
+            year=year, dekad=dekad, local=False
+        )
+        self._download(
+            filepath=filepath, url_filename=url_filename, clobber=clobber
+        )
+
+    def _download(self, filepath: Path, url_filename: str, clobber: bool):
+        local_filename = filepath.stem
+
+        url = self._get_url(filename=url_filename)
+        year, dekad = self._fp_year_dekad(filepath)
+
+        try:
+            resp = urlopen(url)
+        except HTTPError:
+            logger.error(
+                f"No NDVI data available for "
+                f"dekad {dekad} of {year}, skipping."
+            )
+            return
+
+        if filepath.exists():
+            file_time = datetime.fromtimestamp(filepath.stat().st_mtime)
+            url_time_str = resp.headers["last-modified"]
+            url_time = datetime.strptime(
+                url_time_str,
+                # Fri, 27 Mar 2015 08:05:42 GMT
+                "%a, %d %b %Y %X %Z",
+            )
+            if not clobber and url_time <= file_time:
+                logger.info(
+                    f"File {filepath} exists, clobber set to False, and "
+                    "file not been modified since last download. "
+                    "Using existing files."
+                )
+                return filepath
+
+        logger.info(
+            f"Downloading NDVI data for {year}, dekad {dekad} "
+            f"into {filepath}."
+        )
+
+        # open file within memory
+        zf = ZipFile(BytesIO(resp.read()))
+
+        # extract single .tif file from .zip
+        for file in zf.infolist():
+            if file.filename.endswith(".tif"):
+                # rename the file to standardize to name of zip
+                file.filename = f"{local_filename}.tif"
+                zf.extract(file, self._raw_base_dir)
+
+        resp.close()
+        return filepath
+
+    def _process(
+        self,
+        clobber: bool,
+        gdf: gpd.GeoDataFrame,
+        feature_col: str,
+        dates_to_process: list,
+        stat: str,
+        is_percentile: bool,
+        kwargs,
+    ) -> Path:
+        """Process data for particular statistic."""
+        # get processed path for particular statistic
+        percentile_list: Optional[List[int]]
+        if is_percentile:
+            stats_list = None
+            percentile_list = [int(stat)]
+            stat = f"{stat}quant"
+        else:
+            stats_list = [stat]
+            percentile_list = None
+
+        processed_path = self._get_processed_path(
+            feature_col=feature_col, stat=stat
+        )
+
+        if processed_path.is_file():
+            logger.info(
+                f"Processing data from {self._start_date[0]}, "
+                f"dekad {self._start_date[1]} to {self._end_date[0]} "
+                f"dekad {self._end_date[1]} into {processed_path}."
+            )
+
+            (
+                dates_to_process,
+                df_already_processed,
+            ) = self._determine_process_dates(
+                clobber=clobber,
+                filepath=processed_path,
+                dates_to_process=dates_to_process,
+            )
+
+            if not dates_to_process:
+                logger.info(
+                    (
+                        "No new {stat} data to process between "
+                        f"{self._start_date[0]}, "
+                        f"dekad {self._start_date[1]} "
+                        f"and {self._end_date[0]}, "
+                        f"dekad {self._end_date[1]}, "
+                        "set `clobber = True` to re-process this data."
+                    )
+                )
+                return processed_path
+        else:
+            # no dates already processed
+            df_already_processed = pd.DataFrame()
+
+        # process data for necessary dates
+        data = [df_already_processed]
+        for process_date in dates_to_process:
+            da = self.load_raster(process_date)
+            stats = da.aat.compute_raster_stats(
+                gdf=gdf,
+                feature_col=feature_col,
+                stats_list=stats_list,
+                percentile_list=percentile_list,
+                **kwargs,
+            )
+            data.append(stats)
+
+        # join data together and sort
+        df = pd.concat(data).sort_values(by="date").reset_index(drop=True)
+
+        # saving file
+        self._processed_base_dir.mkdir(parents=True, exist_ok=True)
+        df.to_csv(processed_path, index=False)
+
+        return processed_path
+
+    def _determine_process_dates(
+        self, clobber: bool, filepath: Path, dates_to_process: list
+    ) -> Tuple[list, pd.DataFrame]:
+        """Determine dates to process.
+
+        Parameters
+        ----------
+        clobber : bool
+            If True, overwrites existing file
+        filepath : Path
+            Filepath to the existing processed file.
+        dates_to_process : list
+            List of dates to process
+
+        Returns
+        -------
+        Tuple[list, pd.DataFrame]
+            Returns a list of dates to process, filtered
+            based on clobber, and a data frame of existing
+            data to build upon in processing
+
+        Raises
+        ------
+        ValueError
+            Raised if not `clobber` but the statistics
+            and `feature_col` for processing do not
+            match the existing file.
+        """
+        df = self._load(filepath=filepath)
+
+        # get dates that have already been processed
+        dates_already_processed = df[
+            ["year", "dekad", "modified"]
+        ].values.tolist()
+        dates_already_processed = {
+            tuple(d[0:2]): d[2] for d in dates_already_processed
+        }
+
+        if clobber:
+            # remove processed dates from file
+            # so they can be reprocessed
+            keep_rows = [
+                d not in dates_to_process
+                for d in list(dates_already_processed.keys())
+            ]
+            df = df.loc[keep_rows]
+        else:
+            # remove processed dates from dates to process
+            dates_to_process = [
+                d
+                for d in dates_to_process
+                if d not in list(dates_already_processed.keys())
+                or self._get_modified_time(year=d[0], dekad=d[1])
+                > dates_already_processed[d]
+            ]
+
+            # incase dates still processed based on
+            # modified time remove from the df
+            keep_rows = [
+                d not in dates_to_process
+                for d in list(dates_already_processed.keys())
+            ]
+            df = df.loc[keep_rows]
+        return (dates_to_process, df)
+
+    def _get_raw_filename(self, year: int, dekad: int, local: bool) -> str:
+        """Get raw filename (excluding file type suffix).
+
+        Parameters
+        ----------
+        year : int
+            4-digit year
+        dekad : int
+            Dekad
+        local : bool
+            If True, returns filepath for local storage,
+            which includes full 4-digit year and _
+            separating with dekad. If False, filepath
+            corresponds to the zip file stored in the
+            USGS server.
+
+        Returns
+        -------
+        str
+            File path prefix for .zip file at URL and
+            for .tif files stored within the .zip
+        """
+        if local:
+            file_year = f"{year:04}_"
+        else:
+            file_year = f"{year-2000:02}"
+        file_name = (
+            f"{self._datasource_config.area_prefix}{file_year}"
+            f"{dekad:02}{self._data_variable_suffix}"
+        )
+        return file_name
+
+    @staticmethod
+    def _load(filepath: Path, drop_modified: bool = False):
+        """Load processed data.
+
+        Parameters
+        ----------
+        filepath : Path
+            Filepath to processed data
+        drop_modified : bool, default = False
+            Drop modified column, used for merging data frames
+            together.
+
+        Returns
+        -------
+        pd.DataFrame
+            Processed data frame
+        """
+        df = pd.read_csv(filepath, parse_dates=["date", "modified"])
+        if drop_modified:
+            df.drop(["modified"], axis=1, inplace=True)
+        return df
+
+    def _get_raw_path(self, year: int, dekad: int, local: bool) -> Path:
+        """Get raw filepath.
+
+        Parameters
+        ----------
+        year : int
+            4-digit year
+        dekad : int
+            Dekad
+        local : bool
+            If True, returns filepath for local storage,
+            which includes full 4-digit year and _
+            separating with dekad. If False, filepath
+            corresponds to the zip file stored in the
+            USGS server.
+
+        Returns
+        -------
+        Path
+            Path to raw file
+        """
+        filename = self._get_raw_filename(year=year, dekad=dekad, local=local)
+        return self._raw_base_dir / f"{filename}.tif"
+
+    def _get_modified_time(self, year: int, dekad: int) -> datetime:
+        """Get modified time of raw file.
+
+        Used to determine when to re-process
+        or re-download raw raster files.
+
+        Parameters
+        ----------
+        year : int
+            4-digit year
+        dekad : int
+            Dekad
+
+        Returns
+        -------
+        datetime
+            Timestamp of when file was modified.
+        """
+        filepath = self._get_raw_path(year=year, dekad=dekad, local=True)
+        return datetime.fromtimestamp(filepath.stat().st_mtime)
+
+    def _get_processed_filename(self, feature_col: str, stat: str) -> str:
+        """Return processed filename.
+
+        Returns the processed filename. The suffix
+        of the filename is always the ``feature_col``
+        the statistics are aggregated to.
+
+        Returns
+        -------
+        str
+            Processed filename
+        """
+        base_file_name = self._get_processed_base_filename(
+            feature_col=feature_col
+        )
+        return f"{base_file_name}_{stat}.csv"
+
+    def _get_processed_base_filename(self, feature_col: str) -> str:
+        base_file_name = (
+            f"{self._country_config.iso3}"
+            f"_usgs_ndvi_{self._data_variable}"
+            f"_{feature_col}"
+        )
+        return base_file_name
+
+    def _get_processed_path(self, feature_col: str, stat: str) -> Path:
+        return self._processed_base_dir / self._get_processed_filename(
+            feature_col=feature_col, stat=stat
+        )
+
+    def _find_processed_files(self, feature_col: str) -> Generator:
+        base_file_name = self._get_processed_base_filename(
+            feature_col=feature_col
+        )
+        return self._processed_base_dir.glob(f"{base_file_name}_*.csv")
+
+    def _get_url(self, filename) -> str:
+        """Get USGS NDVI URL.
+
+        Parameters
+        ----------
+        filename : str
+            File name string generated for specific year, dekad, and data type
+
+        Returns
+        -------
+        str
+            Download URL string
+        """
+        return (
+            f"https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/fews/web/"
+            f"{self._datasource_config.area_url}/dekadal/emodis"
+            f"/ndvi_c6/{self._data_variable_url}/"
+            f"downloads/dekadal/{filename}.zip"
+        )
+
+    # TODO: potentially move from static method to
+    # wider USGS function repository
+    @staticmethod
+    def _fp_year_dekad(path: Path) -> List[int]:
+        """Extract year and dekad from filepath.
+
+        Parameters
+        ----------
+        path : Path
+            Filepath
+
+        Returns
+        -------
+        list
+            List of year and dekad
+        """
+        filename = path.stem
+        # find two groups, first for year second for dekad
+        regex = re.compile(r"(\d{4})_(\d{2})")
+        return [int(x) for x in regex.findall(filename)[0]]
diff --git a/src/aatoolbox/datasources/usgs/ndvi_products.py b/src/aatoolbox/datasources/usgs/ndvi_products.py
new file mode 100644
index 00000000..18d2f221
--- /dev/null
+++ b/src/aatoolbox/datasources/usgs/ndvi_products.py
@@ -0,0 +1,296 @@
+"""Classes to download and process USGS FEWS NET NDVI data.
+
+Download, process, and load NDVI data published
+in the `USGS FEWS NET data portal
+<https://earlywarning.usgs.gov/fews>`_. Classes available
+to process temporally smoothed NDVI values, percent of
+median values, difference to median values, and difference
+from current value to the previous year's value.
+"""
+
+# TODO: add progress bar
+import logging
+from datetime import date
+from typing import Tuple, Union
+
+from aatoolbox.config.countryconfig import CountryConfig
+from aatoolbox.datasources.usgs.ndvi_base import _UsgsNdvi
+
+logger = logging.getLogger(__name__)
+
+
+class UsgsNdviSmoothed(_UsgsNdvi):
+    """Base class to retrieve smoothed NDVI data.
+
+    The retrieved data is the smoothed NDVI values
+    processed by the USGS. Temporal smoothing is done
+    to adjust for cloud cover and other errors.
+    Data for the 3 most recent dekads is not fully
+    smoothed, and are re-smoothed at the end of the
+    3 dekad period.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    start_date : _DATE_TYPE, default = None
+        Start date. Can be passed as a ``datetime.date``
+        object or a data string in ISO8601 format, and
+        the relevant dekad will be determined. Or pass
+        directly as year-dekad tuple, e.g. (2020, 1).
+        If ``None``, ``start_date`` is set to earliest
+        date with data: 2002, dekad 19.
+    end_date : _DATE_TYPE, default = None
+        End date. Can be passed as a ``datetime.date``
+        object and the relevant dekad will be determined,
+        as a date string in ISO8601 format, or as a
+        year-dekad tuple, i.e. (2020, 1). If ``None``,
+        ``end_date`` is set to ``date.today()``.
+
+    Examples
+    --------
+    >>> from aatoolbox import create_country_config, \
+    ...  CodAB, UsgsNdviSmoothed
+    >>>
+    >>> # Retrieve admin 2 boundaries for Burkina Faso
+    >>> country_config = create_country_config(iso3="bfa")
+    >>> codab = CodAB(country_config=country_config)
+    >>> bfa_admin2 = codab.load(admin_level=2)
+    >>>
+    >>> # setup NDVI
+    >>> bfa_ndvi = UsgsNdviSmoothed(
+    ...     country_config=country_config,
+    ...     start_date=[2020, 1],
+    ...     end_date=[2020, 3]
+    ... )
+    >>> bfa_ndvi.download()
+    >>> bfa_ndvi.process(
+    ...     gdf=bfa_admin2,
+    ...     feature_col="ADM2_FR"
+    ... )
+    >>>
+    >>> # load in processed data
+    >>> df = bfa_ndvi.load(feature_col="ADM2_FR")
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        start_date: Union[date, str, Tuple[int, int], None] = None,
+        end_date: Union[date, str, Tuple[int, int], None] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            data_variable="smoothed",
+            data_variable_suffix="",
+            data_variable_url="temporallysmoothedndvi",
+            start_date=start_date,
+            end_date=end_date,
+        )
+
+
+class UsgsNdviPctMedian(_UsgsNdvi):
+    """Base class to retrieve % of median NDVI.
+
+    The retrieved data is the percent of median NDVI
+    values calculated from 2003 - 2017, as
+    processed by the USGS.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    start_date : _DATE_TYPE, default = None
+        Start date. Can be passed as a ``datetime.date``
+        object or a data string in ISO8601 format, and
+        the relevant dekad will be determined. Or pass
+        directly as year-dekad tuple, e.g. (2020, 1).
+        If ``None``, ``start_date`` is set to earliest
+        date with data: 2002, dekad 19.
+    end_date : _DATE_TYPE, default = None
+        End date. Can be passed as a ``datetime.date``
+        object and the relevant dekad will be determined,
+        as a date string in ISO8601 format, or as a
+        year-dekad tuple, i.e. (2020, 1). If ``None``,
+        ``end_date`` is set to ``date.today()``.
+
+    Examples
+    --------
+    >>> from aatoolbox import create_country_config, \
+    ...  CodAB, UsgsNdviPctMedian
+    >>>
+    >>> # Retrieve admin 2 boundaries for Burkina Faso
+    >>> country_config = create_country_config(iso3="bfa")
+    >>> codab = CodAB(country_config=country_config)
+    >>> bfa_admin2 = codab.load(admin_level=2)
+    >>>
+    >>> # setup NDVI
+    >>> bfa_ndvi = UsgsNdviPctMedian(
+    ...     country_config=country_config,
+    ...     start_date=[2020, 1],
+    ...     end_date=[2020, 3]
+    ... )
+    >>> bfa_ndvi.download()
+    >>> bfa_ndvi.process(
+    ...     gdf=bfa_admin2,
+    ...     feature_col="ADM2_FR"
+    ... )
+    >>>
+    >>> # load in processed data
+    >>> df = bfa_ndvi.load(feature_col="ADM2_FR")
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        start_date: Union[date, str, Tuple[int, int], None] = None,
+        end_date: Union[date, str, Tuple[int, int], None] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            data_variable="percent_median",
+            data_variable_suffix="pct",
+            data_variable_url="percentofmedian",
+            start_date=start_date,
+            end_date=end_date,
+        )
+
+
+class UsgsNdviMedianAnomaly(_UsgsNdvi):
+    """Base class to retrieve NDVI anomaly data.
+
+    The retrieved data is NDVI anomaly data calculated
+    as a subtraction of the median value, based on data
+    from 2003 - 2017, from the
+    current value. Negative values indicate less
+    vegetation than the median, positive values indicate
+    more vegetation.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    start_date : _DATE_TYPE, default = None
+        Start date. Can be passed as a ``datetime.date``
+        object or a data string in ISO8601 format, and
+        the relevant dekad will be determined. Or pass
+        directly as year-dekad tuple, e.g. (2020, 1).
+        If ``None``, ``start_date`` is set to earliest
+        date with data: 2002, dekad 19.
+    end_date : _DATE_TYPE, default = None
+        End date. Can be passed as a ``datetime.date``
+        object and the relevant dekad will be determined,
+        as a date string in ISO8601 format, or as a
+        year-dekad tuple, i.e. (2020, 1). If ``None``,
+        ``end_date`` is set to ``date.today()``.
+
+    Examples
+    --------
+    >>> from aatoolbox import create_country_config, \
+    ...  CodAB, UsgsNdviMedianAnomaly
+    >>>
+    >>> # Retrieve admin 2 boundaries for Burkina Faso
+    >>> country_config = create_country_config(iso3="bfa")
+    >>> codab = CodAB(country_config=country_config)
+    >>> bfa_admin2 = codab.load(admin_level=2)
+    >>>
+    >>> # setup NDVI
+    >>> bfa_ndvi = UsgsNdviMedianAnomaly(
+    ...     country_config=country_config,
+    ...     start_date=[2020, 1],
+    ...     end_date=[2020, 3]
+    ... )
+    >>> bfa_ndvi.download()
+    >>> bfa_ndvi.process(
+    ...     gdf=bfa_admin2,
+    ...     feature_col="ADM2_FR"
+    ... )
+    >>>
+    >>> # load in processed data
+    >>> df = bfa_ndvi.load(feature_col="ADM2_FR")
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        start_date: Union[date, str, Tuple[int, int], None] = None,
+        end_date: Union[date, str, Tuple[int, int], None] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            data_variable="median_anomaly",
+            data_variable_suffix="stmdn",
+            data_variable_url="mediananomaly",
+            start_date=start_date,
+            end_date=end_date,
+        )
+
+
+class UsgsNdviYearDifference(_UsgsNdvi):
+    """Base class to retrieve NDVI year difference data.
+
+    The retrieved data is NDVI yearly difference data,
+    calculated as the subtraction of the previous year's
+    NDVI value from the current year's. Negative
+    values indicate the current vegetation is less
+    than the previous year's, positive that there
+    is more vegetation in the current year.
+
+    Parameters
+    ----------
+    country_config : CountryConfig
+        Country configuration
+    start_date : _DATE_TYPE, default = None
+        Start date. Can be passed as a ``datetime.date``
+        object or a data string in ISO8601 format, and
+        the relevant dekad will be determined. Or pass
+        directly as year-dekad tuple, e.g. (2020, 1).
+        If ``None``, ``start_date`` is set to earliest
+        date with data: 2002, dekad 19.
+    end_date : _DATE_TYPE, default = None
+        End date. Can be passed as a ``datetime.date``
+        object and the relevant dekad will be determined,
+        as a date string in ISO8601 format, or as a
+        year-dekad tuple, i.e. (2020, 1). If ``None``,
+        ``end_date`` is set to ``date.today()``.
+
+    Examples
+    --------
+    >>> from aatoolbox import create_country_config, \
+    ...  CodAB, UsgsNdviDifference
+    >>>
+    >>> # Retrieve admin 2 boundaries for Burkina Faso
+    >>> country_config = create_country_config(iso3="bfa")
+    >>> codab = CodAB(country_config=country_config)
+    >>> bfa_admin2 = codab.load(admin_level=2)
+    >>>
+    >>> # setup NDVI
+    >>> bfa_ndvi = UsgsNdviDifference(
+    ...     country_config=country_config,
+    ...     start_date=[2020, 1],
+    ...     end_date=[2020, 3]
+    ... )
+    >>> bfa_ndvi.download()
+    >>> bfa_ndvi.process(
+    ...     gdf=bfa_admin2,
+    ...     feature_col="ADM2_FR"
+    ... )
+    >>>
+    >>> # load in processed data
+    >>> df = bfa_ndvi.load(feature_col="ADM2_FR")
+    """
+
+    def __init__(
+        self,
+        country_config: CountryConfig,
+        start_date: Union[date, str, Tuple[int, int], None] = None,
+        end_date: Union[date, str, Tuple[int, int], None] = None,
+    ):
+        super().__init__(
+            country_config=country_config,
+            data_variable="difference",
+            data_variable_suffix="dif",
+            data_variable_url="differencepreviousyear",
+            start_date=start_date,
+            end_date=end_date,
+        )
diff --git a/src/aatoolbox/utils/check_file_existence.py b/src/aatoolbox/utils/check_file_existence.py
index 250e813c..f2832c19 100644
--- a/src/aatoolbox/utils/check_file_existence.py
+++ b/src/aatoolbox/utils/check_file_existence.py
@@ -58,11 +58,23 @@ def check_file_existence(
                 " decorator to work."
             )
         ) from err
-    if filepath.exists() and not clobber:
-        logger.info(
-            f"File {filepath} exists and clobber set to False, "
-            f"using existing files"
-        )
+
+    # check existence
+    exist_dict = {True: "exists", False: "does not exist"}
+
+    # check filepath exists -> clobber
+    usage_dict = {
+        True: {True: "overwriting existing", False: "using existing"},
+        False: {True: "downloading new", False: "downloading new"},
+    }
+    fp_exists = filepath.exists()
+
+    logger.info(
+        f"File {filepath} {exist_dict[fp_exists]} and clobber "
+        f"set to {clobber}, {usage_dict[fp_exists][clobber]} file."
+    )
+
+    if fp_exists and not clobber:
         return filepath
     else:
         return wrapped(*args, **kwargs)
diff --git a/src/aatoolbox/utils/dates.py b/src/aatoolbox/utils/dates.py
new file mode 100644
index 00000000..f6a3474a
--- /dev/null
+++ b/src/aatoolbox/utils/dates.py
@@ -0,0 +1,200 @@
+"""Functions for dealing with dates."""
+
+import itertools
+from datetime import date
+from typing import List, Tuple, Union, cast
+
+
+def get_date_from_user_input(input_date: Union[date, str]) -> date:
+    """Return date from string or date input.
+
+    Processes input data in either ``datetime.date``
+    format or as an ISO8601 string. Generates
+    error message if different object provided.
+
+    Parameters
+    ----------
+    input_date : Union[date, str]
+        ``datetime.date`` object or ISO8601 string.
+
+    Returns
+    -------
+    date
+        ``datetime.date``
+    """
+    if isinstance(input_date, date):
+        return input_date
+    try:
+        input_date = date.fromisoformat(input_date)
+    except ValueError as err:
+        raise ValueError(
+            "`date` values passed as a string must follow ISO8601 date "
+            "format: YYYY-MM-DD."
+        ) from err
+    except TypeError as err:
+        raise TypeError(
+            "`date` values must be either an ISO8601 string or "
+            "`datetime.date` object."
+        ) from err
+
+    return input_date
+
+
+def get_dekadal_date(
+    input_date: Union[date, str, Tuple[int, int], None],
+    default_date: Union[date, str, Tuple[int, int], None] = None,
+) -> Tuple[int, int]:
+    """Calculate dekadal date from general input.
+
+    Processes input ``input_date`` and returns two
+    values, the year and dekad. Input can be of
+    format ``datetime.date``, an ISO8601 date
+    string, an already calculated ``(year, dekad)``
+    format date, or ``None``. If ``None``,
+    ``default_date`` is returned. ``default_date``
+    can also be passed in the above formats.
+    """
+    if input_date is None and default_date is not None:
+        input_date = default_date
+
+    # convert date to various values
+    if not isinstance(input_date, (str, date)):
+        input_tuple = cast(Tuple[int, int], input_date)
+        if len(input_tuple) == 2:
+            year, dekad = input_tuple
+            # assert year-dekad values appropriate, not too strict
+            if year < 1000 or year > 9999 or dekad < 1 or dekad > 36:
+                raise ValueError(
+                    f"(year, dekad) tuple ({year}, {dekad}) invalid. "
+                    "Year should be a 4-digit year and dekad between "
+                    "1 and 36."
+                )
+
+        else:
+            raise ValueError(
+                (
+                    "`date` values for dekadal data "
+                    "should be passed in as "
+                    "`datetime.date` objects, tuples "
+                    "of `(year, dekad)` format, or "
+                    "ISO8601 date strings."
+                )
+            )
+
+    else:
+        input_as_date = get_date_from_user_input(input_date)
+        year, dekad = date_to_dekad(input_as_date)
+
+    return year, dekad
+
+
+def dekad_to_date(dekad: Tuple[int, int]) -> date:
+    """Compute date from dekad and year.
+
+    Date computed from dekad and year in
+    datetime object, corresponding to
+    first day of the dekad. This
+    is based on the
+    `common dekadal definition
+    <http://iridl.ldeo.columbia.edu/maproom/Food_Security/Locusts/Regional/Dekadal_Rainfall/index.html>`_
+    of the 1st and 2nd dekad of a month
+    being the first 10 day periods, and
+    the 3rd dekad being the remaining
+    days within that month.
+    """
+    year = dekad[0]
+    month = ((dekad[1] - 1) // 3) + 1
+    day = 10 * ((dekad[1] - 1) % 3) + 1
+    return date(year=year, month=month, day=day)
+
+
+def date_to_dekad(date_obj: date) -> Tuple[int, int]:
+    """Compute dekad and year from date.
+
+    Dekad computed from date. This
+    is based on the
+    `common dekadal definition
+    <http://iridl.ldeo.columbia.edu/maproom/Food_Security/Locusts/Regional/Dekadal_Rainfall/index.html>`_
+    of the 1st and 2nd dekad of a month
+    being the first 10 day periods, and
+    the 3rd dekad being the remaining
+    days within that month.
+    """
+    year = date_obj.year
+    dekad = min((date_obj.day - 1) // 10, 2) + ((date_obj.month - 1) * 3) + 1
+    return (year, dekad)
+
+
+def compare_dekads_lt(
+    dekad1: Tuple[int, int], dekad2: Tuple[int, int]
+) -> bool:
+    """Is year1/dekad1 less than year2/dekad2.
+
+    Compare two pairs of years and dekads,
+    that the first pair are less than the
+    second pair.
+    """
+    y1, d1 = dekad1
+    y2, d2 = dekad2
+    return y1 < y2 or ((y1 == y2) and (d1 < d2))
+
+
+def compare_dekads_lte(
+    dekad1: Tuple[int, int], dekad2: Tuple[int, int]
+) -> bool:
+    """Is year1/dekad1 less than or equal to year2/dekad2.
+
+    Compare two pairs of years and dekads,
+    that the first pair are less than or
+    equal to the second pair.
+    """
+    y1, d1 = dekad1
+    y2, d2 = dekad2
+    return y1 < y2 or ((y1 == y2) and (d1 <= d2))
+
+
+def compare_dekads_gt(
+    dekad1: Tuple[int, int], dekad2: Tuple[int, int]
+) -> bool:
+    """Is year1/dekad1 greater than year2/dekad2.
+
+    Compare two pairs of years and dekads,
+    that the first pair are greater than the
+    second pair.
+    """
+    return compare_dekads_lt(dekad1=dekad2, dekad2=dekad1)
+
+
+def compare_dekads_gte(
+    dekad1: Tuple[int, int], dekad2: Tuple[int, int]
+) -> bool:
+    """Is year1/dekad1 greater than or equal to year2/dekad2.
+
+    Compare two pairs of years and dekads,
+    that the first pair are greater than or
+    equal to the second pair.
+    """
+    return compare_dekads_lte(dekad1=dekad2, dekad2=dekad1)
+
+
+def expand_dekads(
+    dekad1: Tuple[int, int], dekad2: Tuple[int, int]
+) -> List[Tuple[int, int]]:
+    """Expand for all years/dekads between two dates.
+
+    Takes input year and dekads and returns a list
+    of year/dekad lists.
+    """
+    if compare_dekads_gt(dekad1, dekad2):
+        raise ValueError("`dekad1` must be less than or equal to `dekad2`.")
+
+    y1, d1 = dekad1
+    y2, d2 = dekad2
+    year_range = range(y1, y2 + 1)
+    dekad_range = range(1, 37)
+    date_combos = itertools.product(*[year_range, dekad_range])
+
+    def valid(y, d):
+        return not ((y == y1 and d < d1) or (y == y2 and d > d2))
+
+    return [(y, d) for y, d in date_combos if valid(y, d)]
diff --git a/src/aatoolbox/utils/io.py b/src/aatoolbox/utils/io.py
index bc263318..d1f0946d 100644
--- a/src/aatoolbox/utils/io.py
+++ b/src/aatoolbox/utils/io.py
@@ -29,7 +29,9 @@ def download_url(
     """
     save_path.parent.mkdir(exist_ok=True, parents=True)
     # Remove file if already exists
-    save_path.unlink(missing_ok=True)
+    # TODO: Use missing_ok=True once 3.7 is dropped
+    if save_path.exists():
+        save_path.unlink()
 
     # use a session and chunk_size to prevent
     # crashing when downloading large files while
diff --git a/src/aatoolbox/utils/raster.py b/src/aatoolbox/utils/raster.py
index d0c546a1..d5cbe986 100644
--- a/src/aatoolbox/utils/raster.py
+++ b/src/aatoolbox/utils/raster.py
@@ -1,20 +1,19 @@
 """
 Utilities to manipulate and analyze raster data.
 
-This module extends xarray and rioxarray to provide
-additional functionality for raster processing and
-post-processing. The extension is based on the
-guidance for how to `extend xarray
-<http://xarray.pydata.org/en/stable/internals/extending-xarray.html>`_:
+The raster module provides accessor utilities for xarray
+data arrays and datasets accessible using the `aat` accessor.
+These functions are available just by importing directly
+the library using ``import aatoolbox``.
 
-However, since rioxarray already extends xarray, this
+Since rioxarray already extends xarray, this
 modules extensions inherit from the RasterArray and
 RasterDataset extensions respectively. This ensures
 cleaner code in the module as ``rio`` methods are
 available immediately, but also means a couple of
 design decisions are followed.
 
-The xarray.DataArray and xarray.Dataset
+The ``xarray.DataArray`` and ``xarray.Dataset``
 extensions here inherit from rioxarray base classes.
 Thus, methods that are identical for both objects
 are defined in a mixin class ``AatRasterMixin`` which
@@ -22,7 +21,7 @@
 """
 
 import logging
-from typing import Any, Callable, Dict, List, Optional, Union
+from typing import Any, Dict, List, Optional, Union
 
 import geopandas as gpd
 import numpy as np
@@ -32,7 +31,7 @@
 from rioxarray.exceptions import DimensionError, MissingCRS, NoDataInBounds
 from rioxarray.raster_array import RasterArray
 from rioxarray.raster_dataset import RasterDataset
-from rioxarray.rioxarray import _get_data_var_message
+from rioxarray.rioxarray import CRS, _get_data_var_message
 
 logger = logging.getLogger(__name__)
 
@@ -43,9 +42,12 @@ class AatRasterMixin:
     # setting attributes to avoid mypy error, from SO
     # https://stackoverflow.com/questions/53120262/mypy-how-to-
     # ignore-missing-attribute-errors-in-mixins/53228204#53228204
-    _get_obj: Callable
     x_dim: str
     y_dim: str
+    _obj: Union[xr.DataArray, xr.Dataset]
+    _height: int
+    _width: int
+    _crs: CRS
 
     def __init__(self, xarray_obj):
         super().__init__(xarray_obj)
@@ -92,7 +94,7 @@ def longitude_range(self):
         """
         if self._longitude_range is None:
             if self._x_dim is not None:
-                data_obj = self._get_obj(inplace=False)
+                data_obj = self._get_obj_aat(inplace=False)
                 lon_max = data_obj.indexes[self.x_dim].max()
                 if lon_max > 180:
                     self._longitude_range = "360"
@@ -138,7 +140,7 @@ def set_time_dim(
         >>> da.aat.t_dim
         'F'
         """
-        data_obj = self._get_obj(inplace=inplace)
+        data_obj = self._get_obj_aat(inplace=inplace)
         if t_dim not in data_obj.dims:
             raise DimensionError(
                 "Time dimension ({t_dim}) not found."
@@ -193,7 +195,7 @@ def correct_calendar(
         >>> da_crct["t"].attrs["calendar"]
         '360_day'
         """
-        data_obj = self._get_obj(inplace=inplace)
+        data_obj = self._get_obj_aat(inplace=inplace)
         if (
             "calendar" in data_obj[self.t_dim].attrs.keys()
             and data_obj[self.t_dim].attrs["calendar"] == "360"
@@ -216,7 +218,9 @@ def correct_calendar(
 
         return data_obj if not inplace else None
 
-    def invert_coordinates(self, inplace: bool = False) -> xr.DataArray:
+    def invert_coordinates(
+        self, inplace: bool = False
+    ) -> Optional[Union[xr.DataArray, xr.Dataset]]:
         """
         Invert latitude and longitude in data array.
 
@@ -252,28 +256,22 @@ def invert_coordinates(self, inplace: bool = False) -> xr.DataArray:
         ...  coords={"lat":numpy.array([87, 88, 89, 90]),
         ...          "lon":numpy.array([70, 69, 68, 67])}
         ... )
-        >>> da_inv = da.aat.invert_coordinates()
-        >>> da_inv.get_index("lon")
+        >>> da.aat.invert_coordinates(inplace=True)
+        >>> da.get_index("lon")
         Int64Index([67, 68, 69, 70], dtype='int64', name='lon')
-        >>> da_inv.get_index("lat")
+        >>> da.get_index("lat")
         Int64Index([90, 89, 88, 87], dtype='int64', name='lat')
         """
-        data_obj = self._get_obj(inplace=inplace)
+        data_obj = self._get_obj_aat(inplace=inplace)
         lon_inv, lat_inv = self._check_coords_inverted()
 
-        # Flip the raster as necessary (based on the flags)
-        inv_dict = {}
-
         if lon_inv:
             logger.info("Longitude was inverted, reversing coordinates.")
-            inv_dict[self.x_dim] = data_obj[self.x_dim][::-1]
+            data_obj[self.x_dim] = data_obj[self.x_dim][::-1]
 
         if lat_inv:
             logger.info("Latitude was inverted, reversing coordinates.")
-            inv_dict[self.y_dim] = data_obj[self.y_dim][::-1]
-
-        if inv_dict:
-            data_obj = data_obj.reindex(inv_dict)
+            data_obj[self.y_dim] = data_obj[self.y_dim][::-1]
 
         return data_obj if not inplace else None
 
@@ -293,7 +291,7 @@ def _check_coords_inverted(self):
         >>> da.aat._check_coords_inverted()
         (True, False)
         """
-        data_obj = self._get_obj(inplace=False)
+        data_obj = self._get_obj_aat(inplace=False)
         lat = data_obj.get_index(self.y_dim)
         lon = data_obj.get_index(self.x_dim)
         return lon[0] > lon[-1], lat[0] < lat[-1]
@@ -304,7 +302,9 @@ def change_longitude_range(
         """Convert longitude range between -180 to 180 and 0 to 360.
 
         The standard longitude range is from -180 to 180, while some
-        applications use 0 to 360.
+        applications use 0 to 360. This includes
+        ```rasterstats.zonal_stats`` <https://pypi.org/project/rasterstats/>`_,
+        which assumes ranges from 0 to 360.
 
         ``change_longitude_range()`` will convert between the
         two coordinate ranges based on its current state.
@@ -330,8 +330,9 @@ def change_longitude_range(
         --------
         >>> import xarray
         >>> import numpy
+        >>> import pandas
         >>> temp = 15 + 8 * numpy.random.randn(4, 4, 3)
-        >>> precip = 10 * np.random.rand(4, 4, 3)
+        >>> precip = 10 * numpy.random.rand(4, 4, 3)
         >>> ds = xarray.Dataset(
         ...   {
         ...     "temperature": (["lat", "lon", "time"], temp),
@@ -340,18 +341,18 @@ def change_longitude_range(
         ...   coords={
         ...     "lat":numpy.array([87, 88, 89, 90]),
         ...     "lon":numpy.array([5, 120, 199, 360]),
-        ...     "time": pd.date_range("2014-09-06", periods=3)
+        ...     "time": pandas.date_range("2014-09-06", periods=3)
         ...   }
         ... )
         >>> ds_inv = ds.aat.change_longitude_range()
         >>> ds_inv.get_index("lon")
         Int64Index([-161, 0, 5, 120], dtype='int64', name='lon')
         >>> # invert coordinates back to original, in place
-        >>> ds_inv.aat.change_longitude_range(inplace=True)
+        >>> ds_inv.aat.change_longitude_range(to_180_range=False, inplace=True)
         >>> ds_inv.get_index("lon")
         Int64Index([0, 5, 120, 199], dtype='int64', name='lon')
         """
-        data_obj = self._get_obj(inplace=inplace)
+        data_obj = self._get_obj_aat(inplace=inplace)
 
         if to_180_range and self.longitude_range == "360":
             logger.info("Converting longitude from 0 360 to -180 to 180.")
@@ -359,19 +360,54 @@ def change_longitude_range(
             data_obj[self.x_dim] = np.sort(
                 ((data_obj[self.x_dim] + 180) % 360) - 180
             )
-            self._longitude_range = "180"
+            data_obj.aat._longitude_range = "180"
 
         elif not to_180_range and self.longitude_range == "180":
             logger.info("Converting longitude from -180 to 180 to 0 to 360.")
 
             data_obj[self.x_dim] = np.sort(data_obj[self.x_dim] % 360)
-            self._longitude_range = "360"
+            data_obj.aat._longitude_range = "360"
 
         else:
             logger.info("Coordinates already in required range.")
 
         return data_obj if not inplace else None
 
+    def _get_obj_aat(self, inplace: bool) -> Union[xr.DataArray, xr.Dataset]:
+        """
+        Get object to modify.
+
+        Originally directly used the ``rioxarray``
+        _get_obj() method, however, this was silently failing
+        so implementing a similar function directly within
+        this module.
+
+        Parameters
+        ----------
+        inplace : bool
+            If True, returns self
+
+        Returns
+        -------
+        Union[xr.DataArray, xr.Dataset]
+            Object to modify.
+        """
+        if inplace:
+            return self._obj
+
+        obj_copy = self._obj.copy(deep=True)
+
+        # preserve attribute information
+        obj_copy.rio._x_dim = self._x_dim
+        obj_copy.rio._y_dim = self._y_dim
+        obj_copy.rio._width = self._width
+        obj_copy.rio._height = self._height
+        obj_copy.rio._crs = self._crs
+        obj_copy.aat._t_dim = self._t_dim
+        obj_copy.aat._longitude_range = self._longitude_range
+
+        return obj_copy
+
 
 @xr.register_dataarray_accessor("aat")
 class AatRasterArray(AatRasterMixin, RasterArray):
@@ -449,10 +485,10 @@ def compute_raster_stats(
         0       3.0  1.5811388300841898        1        5     12.0          4  area_a  # noqa: E501
         1       4.5                 1.5        3        6      9.0          2  area_b  # noqa: E501
         """
-        data_obj = self._get_obj(inplace=False)
-        if data_obj.rio.crs is None:
+        if self._obj.rio.crs is None:
             raise MissingCRS("No CRS found, set CRS before computation.")
 
+        data_obj = self._get_obj_aat(inplace=False)
         df_list = []
 
         if stats_list is None:
@@ -491,14 +527,14 @@ def compute_raster_stats(
                     kwargs["min_count"] = 1
                 grid_stat = getattr(da_clip, stat)(
                     dim=[self.x_dim, self.y_dim], **kwargs
-                ).rename(f"{stat}_{feature_col}")
+                ).rename(stat)
                 grid_stat_all.append(grid_stat)
 
             if percentile_list is not None:
                 grid_quant = [
                     da_clip.quantile(quant / 100, dim=[self.x_dim, self.y_dim])
                     .drop("quantile")
-                    .rename(f"{quant}quant_{feature_col}")
+                    .rename(f"{quant}quant")
                     for quant in percentile_list
                 ]
                 grid_stat_all.extend(grid_quant)
@@ -656,7 +692,7 @@ def compute_raster_stats(
         ...    gdf=gdf,
         ...    feature_col="name"
         ... ) # doctest: +SKIP
-        [  mean_name            std_name min_name max_name sum_name count_name    name # noqa: E501
+        [       mean                 std      min      max      sum      count    name # noqa: E501
         0       3.0  1.5811388300841898        1        5     12.0          4  area_a  # noqa: E501
         1       4.5                 1.5        3        6      9.0          2  area_b] # noqa: E501
         """
diff --git a/tox.ini b/tox.ini
index 6dd6ddc1..7b2a7cb2 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,11 +1,11 @@
 [tox]
 isolated_build = true
 envlist =
-    py36
     py37
     py38
     py39
     py310
+    py311
 
 [testenv]
 wheel = true
@@ -21,8 +21,8 @@ commands =
 
 [gh-actions]
 python =
-    3.6: py36
     3.7: py37
     3.8: py38
     3.9: py39
     3.10: py310
+    3.11: py311


Test patch:
diff --git a/.github/workflows/run-python-tests.yaml b/.github/workflows/run-python-tests.yaml
index 4f2cdebd..61c56eca 100644
--- a/.github/workflows/run-python-tests.yaml
+++ b/.github/workflows/run-python-tests.yaml
@@ -1,5 +1,9 @@
 name: tests
-on: push
+on:
+  schedule:
+    - cron: 0 0 * * *
+  push:
+
 jobs:
   build-and-test:
     name: Run unit tests
@@ -7,7 +11,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ['3.6', '3.7', '3.8', '3.9', '3.10']
+        python-version: ['3.7', '3.8', '3.9', '3.10', '3.11']
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}
@@ -34,7 +38,7 @@ jobs:
         path: .tox/test-results.xml
     - name: Upload coverage results to Codecov
       uses: codecov/codecov-action@v2
-      if: success() && matrix.python-version == 3.10
+      if: success() && matrix.python-version == 3.11
 
   publish-test-results:
     name: Publish unit test results
@@ -47,7 +51,7 @@ jobs:
         with:
           path: artifacts
       - name: Publish test results to GitHub
-        uses: EnricoMi/publish-unit-test-result-action@v1
+        uses: EnricoMi/publish-unit-test-result-action@v2
         with:
           github_token: ${{ secrets.GITHUB_TOKEN }}
           files: artifacts/**/*.xml
diff --git a/requirements/requirements-tests.txt b/requirements/requirements-tests.txt
index 90ea17fc..790a66ed 100644
--- a/requirements/requirements-tests.txt
+++ b/requirements/requirements-tests.txt
@@ -1,22 +1,28 @@
-attrs==21.4.0
+#
+# This file is autogenerated by pip-compile with python 3.11
+# To update, run:
+#
+#    pip-compile --output-file=requirements/requirements-tests.txt requirements/tests.in
+#
+attrs==22.1.0
     # via pytest
-cftime==1.6.0
+cftime==1.6.2
     # via
     #   -r requirements/tests.in
     #   netcdf4
-coverage[toml]==6.3.2
+coverage[toml]==6.5.0
     # via pytest-cov
-distlib==0.3.4
+distlib==0.3.6
     # via virtualenv
-filelock==3.6.0
+filelock==3.8.0
     # via
     #   tox
     #   virtualenv
 iniconfig==1.1.1
     # via pytest
-netcdf4==1.5.8
+netcdf4==1.6.2
     # via -r requirements/tests.in
-numpy==1.22.3
+numpy==1.23.5
     # via
     #   cftime
     #   netcdf4
@@ -24,38 +30,28 @@ packaging==21.3
     # via
     #   pytest
     #   tox
-platformdirs==2.5.2
+platformdirs==2.5.4
     # via virtualenv
 pluggy==1.0.0
     # via
     #   pytest
     #   tox
 py==1.11.0
-    # via
-    #   pytest
-    #   tox
-pyparsing==3.0.8
+    # via tox
+pyparsing==3.0.9
     # via packaging
-pytest==7.1.1
+pytest==7.2.0
     # via
     #   -r requirements/tests.in
     #   pytest-cov
     #   pytest-mock
-pytest-cov==3.0.0
+pytest-cov==4.0.0
     # via -r requirements/tests.in
-pytest-mock==3.7.0
+pytest-mock==3.10.0
     # via -r requirements/tests.in
 six==1.16.0
-    # via
-    #   tox
-    #   virtualenv
-toml==0.10.2
     # via tox
-tomli==2.0.1
-    # via
-    #   coverage
-    #   pytest
-tox==3.25.0
+tox==3.27.1
     # via -r requirements/tests.in
-virtualenv==20.14.1
+virtualenv==20.16.7
     # via tox
diff --git a/tests/datasources/conftest.py b/tests/datasources/conftest.py
index af3239c6..945dec6e 100644
--- a/tests/datasources/conftest.py
+++ b/tests/datasources/conftest.py
@@ -1,7 +1,7 @@
 """Fixtures for all pipeline-related tests."""
 import pytest
 
-from aatoolbox import create_custom_country_config
+from aatoolbox import GeoBoundingBox, create_custom_country_config
 
 CONFIG_FILE = "tests/datasources/fake_config.yaml"
 ISO3 = "abc"
@@ -27,6 +27,13 @@ def mock_country_config():
     return create_custom_country_config(filepath=CONFIG_FILE)
 
 
+@pytest.fixture
+def geo_bounding_box():
+    """Input GeoBoundingBox to use."""
+    gbb = GeoBoundingBox(lat_max=1.0, lat_min=-2.2, lon_max=3.3, lon_min=-4.4)
+    return gbb
+
+
 def pytest_configure(config):
     """Create custom markers to add to tests."""
     config.addinivalue_line(
diff --git a/tests/datasources/fake_config.yaml b/tests/datasources/fake_config.yaml
index 707d14e9..d6ea950b 100644
--- a/tests/datasources/fake_config.yaml
+++ b/tests/datasources/fake_config.yaml
@@ -3,5 +3,16 @@ codab:
   hdx_dataset_name: fake_hdx_dataset_name
   layer_base_name: fake_layer_base_name_level{admin_level}
   admin_level_max: 2
+  admin2_name: admin2_custom_name
 fewsnet:
   region_name: east-africa
+glofas:
+  reporting_points:
+    - name: Station1
+      lon: 2.5
+      lat: 0.5
+    - name: Station2
+      lon: -3.5
+      lat: -1.5
+usgs_ndvi:
+  area_name: east-africa
diff --git a/tests/datasources/test_chirps.py b/tests/datasources/test_chirps.py
new file mode 100644
index 00000000..ca49fab0
--- /dev/null
+++ b/tests/datasources/test_chirps.py
@@ -0,0 +1,480 @@
+"""Tests for the Chirps module."""
+import calendar
+from datetime import date
+from pathlib import Path
+
+import cftime
+import numpy as np
+import pytest
+import xarray as xr
+from xarray.coding.cftimeindex import CFTimeIndex
+
+from aatoolbox import ChirpsDaily, ChirpsMonthly, GeoBoundingBox
+
+DATASOURCE_BASE_DIR = "chirps"
+
+START_DATE = date(year=2020, month=7, day=15)
+END_DATE = date(year=2020, month=8, day=1)
+CURRENT_DATE = date(year=2022, month=5, day=31)
+PAST_DATE = date(year=1970, month=1, day=1)
+FUTURE_DATE = date(year=2100, month=10, day=2)
+
+START_YEAR = "2020"
+START_MONTH = "07"
+START_MONTH_NAME = "Jul"
+START_DAY = "15"
+
+END_YEAR = "2020"
+END_MONTH = "08"
+END_MONTH_NAME = "Aug"
+END_DAY = "01"
+
+
+@pytest.fixture
+def mock_chirps(mocker, mock_country_config):
+    """Create Chirps class with mock country config."""
+    geo_bounding_box = GeoBoundingBox(
+        lat_max=6, lat_min=3.2, lon_max=2, lon_min=-3
+    )
+
+    mocker.patch(
+        (
+            "aatoolbox.datasources.chirps.chirps.ChirpsDaily"
+            "._get_last_available_date"
+        ),
+        return_value=CURRENT_DATE,
+    )
+
+    mocker.patch(
+        (
+            "aatoolbox.datasources.chirps.chirps.ChirpsMonthly"
+            "._get_last_available_date"
+        ),
+        return_value=CURRENT_DATE,
+    )
+
+    def _mock_chirps(
+        frequency: str = "daily",
+        resolution: float = 0.05,
+        start_date: date = START_DATE,
+        end_date: date = END_DATE,
+    ):
+        if frequency == "daily":
+            chirps = ChirpsDaily(
+                country_config=mock_country_config,
+                geo_bounding_box=geo_bounding_box,
+                resolution=resolution,
+                start_date=start_date,
+                end_date=end_date,
+            )
+        else:
+            chirps = ChirpsMonthly(
+                country_config=mock_country_config,
+                geo_bounding_box=geo_bounding_box,
+                start_date=start_date,
+                end_date=end_date,
+            )
+        return chirps
+
+    return _mock_chirps
+
+
+@pytest.fixture
+def mock_xr_open_multiple_dataset(mocker):
+    """Mock GeoPandas file reading function."""
+    ds = xr.Dataset()
+    return mocker.patch(
+        "aatoolbox.datasources.chirps.chirps.xr.open_mfdataset",
+        return_value=ds,
+    )
+
+
+@pytest.fixture
+def mock_dataset_to_netcdf(mocker):
+    """Mock GeoPandas file writing function."""
+    return mocker.patch(
+        "aatoolbox.datasources.chirps.chirps.xr.Dataset.to_netcdf",
+    )
+
+
+@pytest.fixture
+def mock_download(mocker, mock_chirps):
+    """Create mock for download method."""
+    download_mock = mocker.patch(
+        "aatoolbox.datasources.chirps.chirps._Chirps._download"
+    )
+
+    def _mock_download(
+        frequency: str,
+    ):
+        chirps = mock_chirps(
+            frequency=frequency,
+            start_date=START_DATE,
+            end_date=END_DATE,
+        )
+        chirps.download()
+        args_download = download_mock.call_args_list
+        url_list = [k["url"] for (_, k) in args_download]
+        filepath_list = [k["filepath"] for (_, k) in args_download]
+        return url_list, filepath_list
+
+    return _mock_download
+
+
+def test_valid_arguments_class(mocker, mock_country_config):
+    """Test for resolution and wrong dates in initialisation class."""
+    mocker.patch(
+        (
+            "aatoolbox.datasources.chirps.chirps.ChirpsDaily"
+            "._get_last_available_date"
+        ),
+        return_value=CURRENT_DATE,
+    )
+
+    geo_bounding_box = GeoBoundingBox(
+        lat_max=6, lat_min=3.2, lon_max=2, lon_min=-3
+    )
+    with pytest.raises(ValueError):
+        ChirpsDaily(
+            country_config=mock_country_config,
+            geo_bounding_box=geo_bounding_box,
+            resolution=0.10,
+            start_date=START_DATE,
+            end_date=END_DATE,
+        )
+        ChirpsDaily(
+            country_config=mock_country_config,
+            geo_bounding_box=geo_bounding_box,
+            resolution=0.10,
+            start_date=START_DATE,
+            end_date=FUTURE_DATE,
+        )
+        ChirpsDaily(
+            country_config=mock_country_config,
+            geo_bounding_box=geo_bounding_box,
+            resolution=0.10,
+            start_date=PAST_DATE,
+            end_date=END_DATE,
+        )
+
+
+def test_download_monthly(
+    mock_aa_data_dir, mock_country_config, mock_download
+):
+    """Test of call download for monthly data."""
+    url_list, filepath_list = mock_download(frequency="monthly")
+
+    url_list_control = [
+        (
+            "https://iridl.ldeo.columbia.edu/SOURCES/.UCSB/"
+            ".CHIRPS/.v2p0/.monthly/.global/.precipitation/"
+            f"X/%28-3.0%29%282.0%29RANGEEDGES/"
+            f"Y/%286.0%29%283.2%29RANGEEDGES/"
+            f"T/%28{MONTH_NAME}%20{START_YEAR}%29"
+            f"%28{MONTH_NAME}%20{START_YEAR}"
+            "%29RANGEEDGES/data.nc"
+        )
+        for MONTH_NAME in [
+            calendar.month_abbr[month]
+            for month in range(int(START_MONTH), int(END_MONTH) + 1)
+        ]
+    ]
+
+    filepath_list_control = [
+        mock_aa_data_dir
+        / (
+            f"public/raw/{mock_country_config.iso3}/{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_monthly_{START_YEAR}_{MONTH:02d}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+        for MONTH in range(int(START_MONTH), int(END_MONTH) + 1)
+    ]
+
+    assert url_list == url_list_control
+    assert filepath_list == filepath_list_control
+
+
+def test_download_daily(mock_aa_data_dir, mock_country_config, mock_download):
+    """Test of call download for daily data."""
+    url_list, filepath_list = mock_download(frequency="daily")
+
+    url_list_control = [
+        (
+            "https://iridl.ldeo.columbia.edu/SOURCES/.UCSB/"
+            ".CHIRPS/.v2p0/.daily-improved/.global/.0p05/.prcp/"
+            f"X/%28-3.0%29%282.0%29RANGEEDGES/"
+            f"Y/%286.0%29%283.2%29RANGEEDGES/"
+            f"T/%28{DAY:02d}%20{START_MONTH_NAME}%20{START_YEAR}"
+            f"%29%28{DAY:02d}%20{START_MONTH_NAME}%20{START_YEAR}"
+            "%29RANGEEDGES/data.nc"
+        )
+        for DAY in range(int(START_DAY), 32)
+    ]
+
+    url_list_control.append(
+        "https://iridl.ldeo.columbia.edu/SOURCES/.UCSB/"
+        ".CHIRPS/.v2p0/.daily-improved/.global/.0p05/.prcp/"
+        f"X/%28-3.0%29%282.0%29RANGEEDGES/"
+        f"Y/%286.0%29%283.2%29RANGEEDGES/"
+        f"T/%2801%20{END_MONTH_NAME}%20{START_YEAR}"
+        f"%29%2801%20{END_MONTH_NAME}%20{START_YEAR}"
+        "%29RANGEEDGES/data.nc"
+    )
+
+    filepath_list_control = [
+        mock_aa_data_dir
+        / (
+            f"public/raw/{mock_country_config.iso3}/{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_daily_{START_YEAR}_{START_MONTH}_{DAY}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+        for DAY in range(int(START_DAY), 32)
+    ]
+
+    filepath_list_control.append(
+        mock_aa_data_dir
+        / (
+            f"public/raw/{mock_country_config.iso3}/{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_daily_{START_YEAR}_{END_MONTH}_{END_DAY}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+    )
+
+    assert url_list == url_list_control
+    assert filepath_list == filepath_list_control
+
+
+def test_process_monthly(
+    mocker,
+    mock_chirps,
+    mock_aa_data_dir,
+    mock_country_config,
+    mock_dataset_to_netcdf,
+):
+    """Test process monthl data."""
+    ds = xr.DataArray(
+        np.reshape(a=np.arange(8), newshape=(2, 2, 2)),
+        dims=("X", "Y", "T"),
+        coords={
+            "X": [1, 2],
+            "Y": [2, -3],
+            "T": [685.5, 686.5],
+        },
+    ).to_dataset(name="prcp")
+
+    ds["T"].attrs["calendar"] = "360"
+    ds["T"].attrs["units"] = "months since 1960-01-01"
+
+    mock_xr_open_dataset = mocker.patch(
+        "aatoolbox.datasources.chirps.chirps.xr.open_dataset",
+        return_value=ds,
+    )
+
+    input_filepath_list_control = [
+        mock_aa_data_dir
+        / (
+            f"public/raw/{mock_country_config.iso3}/"
+            f"{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_monthly_{START_YEAR}_{MONTH:02d}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+        for MONTH in range(int(START_MONTH), int(END_MONTH) + 1)
+    ]
+
+    output_filepath_list_control = [
+        Path(str(f).replace("raw", "processed"))
+        for f in input_filepath_list_control
+    ]
+
+    chirps = mock_chirps(frequency="monthly")
+    processed_path = chirps.process()
+
+    args_input_process = mock_xr_open_dataset.call_args_list
+    input_filepath_list = [k[0] for (k, _) in args_input_process]
+
+    args_output_process = mock_dataset_to_netcdf.call_args_list
+    output_filepath_list = [k["path"] for (_, k) in args_output_process]
+    output_dataset_list = [k[0] for (k, _) in args_output_process]
+    output_ds = output_dataset_list[0]
+
+    assert input_filepath_list == input_filepath_list_control
+    assert output_filepath_list == output_filepath_list_control
+
+    assert processed_path == (
+        mock_aa_data_dir
+        / f"public/processed/{mock_country_config.iso3}/{DATASOURCE_BASE_DIR}/"
+    )
+
+    expected_f = CFTimeIndex(
+        [
+            cftime.Datetime360Day(year=2017, month=2, day=16),
+            cftime.Datetime360Day(year=2017, month=3, day=16),
+        ]
+    )
+
+    assert np.array_equal(output_ds.X.values, [1, 2])
+    assert np.array_equal(output_ds.Y.values, [2, -3])
+    assert output_ds.get_index("T").equals(expected_f)
+    assert np.array_equal(
+        output_ds.precipitation.values,
+        np.reshape(a=np.arange(8), newshape=(2, 2, 2)),
+    )
+
+
+def test_process_daily(
+    mocker,
+    mock_chirps,
+    mock_aa_data_dir,
+    mock_country_config,
+    mock_dataset_to_netcdf,
+):
+    """Test process daily data."""
+    ds = xr.DataArray(
+        np.reshape(a=np.arange(8), newshape=(2, 2, 2)),
+        dims=("X", "Y", "T"),
+        coords={
+            "X": [1, 2],
+            "Y": [2, -3],
+            "T": [2454397.0, 2454398.0],
+        },
+    ).to_dataset(name="prcp")
+
+    ds["T"].attrs["calendar"] = "standard"
+    ds["T"].attrs["units"] = "julian_day"
+
+    mock_xr_open_dataset = mocker.patch(
+        "aatoolbox.datasources.chirps.chirps.xr.open_dataset",
+        return_value=ds,
+    )
+
+    input_filepath_list_control = [
+        mock_aa_data_dir
+        / (
+            f"public/raw/{mock_country_config.iso3}/"
+            f"{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_daily_{START_YEAR}_{START_MONTH}_{DAY}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+        for DAY in range(int(START_DAY), 32)
+    ]
+
+    input_filepath_list_control.append(
+        mock_aa_data_dir
+        / (
+            f"public/raw/{mock_country_config.iso3}/"
+            f"{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_daily_{START_YEAR}_{END_MONTH}_{END_DAY}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+    )
+
+    output_filepath_list_control = [
+        Path(str(f).replace("raw", "processed"))
+        for f in input_filepath_list_control
+    ]
+
+    chirps = mock_chirps(frequency="daily")
+    processed_path = chirps.process()
+
+    args_input_process = mock_xr_open_dataset.call_args_list
+    input_filepath_list = [k[0] for (k, _) in args_input_process]
+
+    args_output_process = mock_dataset_to_netcdf.call_args_list
+    output_filepath_list = [k["path"] for (_, k) in args_output_process]
+    output_dataset_list = [k[0] for (k, _) in args_output_process]
+    output_ds = output_dataset_list[0]
+
+    assert input_filepath_list == input_filepath_list_control
+    assert output_filepath_list == output_filepath_list_control
+
+    assert processed_path == (
+        mock_aa_data_dir
+        / f"public/processed/{mock_country_config.iso3}/{DATASOURCE_BASE_DIR}/"
+    )
+
+    expected_f = CFTimeIndex(
+        [
+            cftime.DatetimeGregorian(year=2007, month=10, day=23, hour=12),
+            cftime.DatetimeGregorian(year=2007, month=10, day=24, hour=12),
+        ]
+    )
+
+    assert np.array_equal(output_ds.X.values, [1, 2])
+    assert np.array_equal(output_ds.Y.values, [2, -3])
+    assert output_ds.get_index("T").equals(expected_f)
+    assert np.array_equal(
+        output_ds.precipitation.values,
+        np.reshape(a=np.arange(8), newshape=(2, 2, 2)),
+    )
+
+
+def test_chirps_load_monthly(
+    mock_xr_open_multiple_dataset,
+    mock_chirps,
+    mock_aa_data_dir,
+    mock_country_config,
+):
+    """Test load monthly data."""
+    filepath_list_control = [
+        mock_aa_data_dir
+        / (
+            f"public/processed/{mock_country_config.iso3}/"
+            f"{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_monthly_{START_YEAR}_{MONTH:02d}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+        for MONTH in range(int(START_MONTH), int(END_MONTH) + 1)
+    ]
+
+    chirps = mock_chirps(frequency="monthly")
+    chirps.load()
+    args_download = mock_xr_open_multiple_dataset.call_args
+    filepath_list = args_download[0][0]
+
+    ds = mock_xr_open_multiple_dataset()
+
+    assert ds.attrs["included_files"] == [
+        filepath.stem for filepath in filepath_list
+    ]
+    assert filepath_list == filepath_list_control
+
+
+def test_chirps_load_daily(
+    mock_xr_open_multiple_dataset,
+    mock_chirps,
+    mock_aa_data_dir,
+    mock_country_config,
+):
+    """Test load daily data."""
+    filepath_list_control = [
+        mock_aa_data_dir
+        / (
+            f"public/processed/{mock_country_config.iso3}/"
+            f"{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_daily_{START_YEAR}_{START_MONTH}_{DAY}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+        for DAY in range(int(START_DAY), 32)
+    ]
+
+    filepath_list_control.append(
+        mock_aa_data_dir
+        / (
+            f"public/processed/{mock_country_config.iso3}/"
+            f"{DATASOURCE_BASE_DIR}/"
+            f"abc_chirps_daily_{START_YEAR}_{END_MONTH}_{END_DAY}_"
+            "r0.05_Np6Sp3Ep2Wm3.nc"
+        )
+    )
+
+    chirps = mock_chirps(frequency="daily")
+    chirps.load()
+    args_download = mock_xr_open_multiple_dataset.call_args
+    filepath_list = args_download[0][0]
+
+    ds = mock_xr_open_multiple_dataset()
+
+    assert ds.attrs["included_files"] == [
+        filepath.stem for filepath in filepath_list
+    ]
+    assert filepath_list == filepath_list_control
diff --git a/tests/datasources/test_codab.py b/tests/datasources/test_codab.py
index 894f2f07..02640fee 100644
--- a/tests/datasources/test_codab.py
+++ b/tests/datasources/test_codab.py
@@ -1,6 +1,4 @@
 """Test COD AB methods."""
-from pathlib import Path
-
 import pytest
 
 from aatoolbox import CodAB
@@ -41,11 +39,20 @@ def test_codab_load_admin_level(
 ):
     """Test that load_codab retrieves expected file and layer name."""
     codab = CodAB(country_config=mock_country_config)
-    admin_level = 2
-    expected_layer_name = mock_country_config.codab.layer_base_name.format(
-        admin_level=admin_level
+
+    # First checking layer_base_name
+    expected_layer_name = "fake_layer_base_name_level1"
+    codab.load(admin_level=1)
+
+    gpd_read_file.assert_called_with(
+        f"zip://{mock_aa_data_dir}/public/raw/{mock_country_config.iso3}/"
+        f"{DATASOURCE_BASE_DIR}/{mock_country_config.iso3}_"
+        f"{DATASOURCE_BASE_DIR}.shp.zip/{expected_layer_name}"
     )
-    codab.load(admin_level=admin_level)
+
+    # Then checking custom name
+    expected_layer_name = "admin2_custom_name"
+    codab.load(admin_level=2)
 
     gpd_read_file.assert_called_with(
         f"zip://{mock_aa_data_dir}/public/raw/{mock_country_config.iso3}/"
@@ -87,8 +94,9 @@ def test_codab_load_fail(mock_country_config):
     """Test raises file not found error when load fails."""
     codab = CodAB(country_config=mock_country_config)
     # Remove file if it exists
+    # TODO: Use missing_ok=True once 3.7 is dropped
     if codab._raw_filepath.exists():
-        Path.unlink(codab._raw_filepath)
+        codab._raw_filepath.unlink()
     with pytest.raises(FileNotFoundError) as excinfo:
         codab.load(admin_level=0)
     assert (
diff --git a/tests/datasources/test_datasource.py b/tests/datasources/test_datasource.py
index ba13c1ff..8c1b609c 100644
--- a/tests/datasources/test_datasource.py
+++ b/tests/datasources/test_datasource.py
@@ -4,7 +4,7 @@
 from aatoolbox.datasources.datasource import DataSource
 
 
-class TestClass(DataSource):
+class DataSourceTesting(DataSource):
     """Class to inherit from DataSource to enable testing of DataSource."""
 
     def __init__(
@@ -40,19 +40,28 @@ def load(self):
 
 def test_config_attribute_name_validator(mock_country_config):
     """Test that correctly checked that datasource_name in the config."""
-    with pytest.raises(AttributeError):
-        TestClass(
-            country_config=mock_country_config,
-            datasource_base_dir="fake_dir_name",
-            config_datasource_name="fake_config_attribute",
-        )
-
-    testclass_valid_attribute_name = TestClass(
+    # Test that attributes are passed to _datasource_config
+    testclass_valid_attribute_name = DataSourceTesting(
         country_config=mock_country_config,
-        datasource_base_dir="fewsnet",
+        datasource_base_dir="fake_dir_name",
         config_datasource_name="fewsnet",
     )
     assert (
         testclass_valid_attribute_name._datasource_config.region_name
         == "east-africa"
     )
+    # Test that non-existant attribute throws error
+    with pytest.raises(AttributeError):
+        DataSourceTesting(
+            country_config=mock_country_config,
+            datasource_base_dir="fake_dir_name",
+            config_datasource_name="fake_config_attribute",
+        )
+    # Test the None attribute throws error
+    mock_country_config.fewsnet = None
+    with pytest.raises(AttributeError):
+        DataSourceTesting(
+            country_config=mock_country_config,
+            datasource_base_dir="fake_dir_name",
+            config_datasource_name="fewsnet",
+        )
diff --git a/tests/datasources/test_fewsnet.py b/tests/datasources/test_fewsnet.py
index ea3effbf..1fb04936 100644
--- a/tests/datasources/test_fewsnet.py
+++ b/tests/datasources/test_fewsnet.py
@@ -37,13 +37,6 @@ def test_no_iso2(mock_country_config):
         FewsNet(country_config=mock_country_config)
 
 
-def test_no_fewsnet_config(mock_country_config):
-    """Test AttributeError when no FEWS NET config."""
-    mock_country_config.fewsnet = None
-    with pytest.raises(AttributeError):
-        FewsNet(country_config=mock_country_config)
-
-
 def test_download_country(
     mock_aa_data_dir, mock_country_config, mock_download_call, mocker
 ):
diff --git a/tests/datasources/test_glofas.py b/tests/datasources/test_glofas.py
new file mode 100644
index 00000000..7058975b
--- /dev/null
+++ b/tests/datasources/test_glofas.py
@@ -0,0 +1,99 @@
+"""Test the GloFAS instantiation."""
+from datetime import date
+
+import pytest
+
+from aatoolbox import GlofasForecast, GlofasReanalysis, GlofasReforecast
+
+
+def test_reanalysis_dates(mock_country_config, geo_bounding_box):
+    """Test date range behaviour of reanalysis."""
+
+    def glofas_reanalysis(end_date: date, start_date: date = None):
+        return GlofasReanalysis(
+            country_config=mock_country_config,
+            geo_bounding_box=geo_bounding_box,
+            start_date=start_date,
+            end_date=end_date,
+        )
+
+    # Try using a string as a date
+    good_end_date = "2022-11-15"
+    # These should not throw an error
+    glofas_reanalysis(end_date=good_end_date)
+    # End date too far in future
+    glofas_future = glofas_reanalysis(end_date=date(year=3000, month=1, day=1))
+    assert glofas_future._end_date == date.today()
+    # Start date too early
+    glofas_past = glofas_reanalysis(
+        start_date=date(year=1800, month=1, day=1), end_date=good_end_date
+    )
+    assert glofas_past._start_date == date(year=1979, month=1, day=1)
+    # End date > start date
+    with pytest.raises(ValueError):
+        glofas_reanalysis(
+            start_date=date(year=2020, month=1, day=2),
+            end_date=date(year=2020, month=1, day=1),
+        )
+
+
+def test_forecast_dates(mock_country_config, geo_bounding_box):
+    """Test date range behaviour of forecast."""
+
+    def glofas_forecast(end_date: date, start_date: date = None):
+        return GlofasForecast(
+            country_config=mock_country_config,
+            geo_bounding_box=geo_bounding_box,
+            start_date=start_date,
+            end_date=end_date,
+            leadtime_max=15,
+        )
+
+    # Try using a string as a date
+    good_end_date = "2022-11-15"
+    # These should not throw an error
+    glofas_forecast(end_date=good_end_date)
+    # End date too far in future
+    glofas_future = glofas_forecast(end_date=date(year=3000, month=1, day=1))
+    assert glofas_future._end_date == date.today()
+    # Start date too early
+    glofas_past = glofas_forecast(
+        start_date=date(year=1800, month=1, day=1), end_date=good_end_date
+    )
+    assert glofas_past._start_date == date(year=2021, month=5, day=26)
+    # End date > start date
+    with pytest.raises(ValueError):
+        glofas_forecast(
+            start_date=date(year=2020, month=1, day=2),
+            end_date=date(year=2020, month=1, day=1),
+        )
+
+
+def test_reforecast_dates(mock_country_config, geo_bounding_box):
+    """Test date range behaviour of reforecast."""
+
+    def glofas_reforecast(start_date: date = None, end_date: date = None):
+        return GlofasReforecast(
+            country_config=mock_country_config,
+            geo_bounding_box=geo_bounding_box,
+            start_date=start_date,
+            end_date=end_date,
+            leadtime_max=15,
+        )
+
+    # These should not throw an error
+    glofas_reforecast()
+    # Try using a string as a date
+    glofas_reforecast(end_date="2010-01-01")
+    # End date too far in future
+    glofas_future = glofas_reforecast(end_date=date(year=3000, month=1, day=1))
+    assert glofas_future._end_date == date(year=2018, month=12, day=31)
+    # Start date too early
+    glofas_past = glofas_reforecast(start_date=date(year=1800, month=1, day=1))
+    assert glofas_past._start_date == date(year=1999, month=1, day=1)
+    # End date > start date
+    with pytest.raises(ValueError):
+        glofas_reforecast(
+            start_date=date(year=2020, month=1, day=2),
+            end_date=date(year=2020, month=1, day=1),
+        )
diff --git a/tests/datasources/test_glofas_download.py b/tests/datasources/test_glofas_download.py
new file mode 100644
index 00000000..88485f39
--- /dev/null
+++ b/tests/datasources/test_glofas_download.py
@@ -0,0 +1,206 @@
+"""Tests for GloFAS data download and processing."""
+from datetime import date
+from pathlib import Path
+from unittest import mock
+
+import pytest
+
+from aatoolbox import (
+    GeoBoundingBox,
+    GlofasForecast,
+    GlofasReanalysis,
+    GlofasReforecast,
+)
+
+
+@pytest.fixture
+def mock_retrieve(mocker):
+    """Mock only the retrieve method of the Client."""
+    mocker.patch(
+        "aatoolbox.datasources.glofas.glofas.cdsapi.Client.__init__",
+        return_value=None,
+    )
+    return mocker.patch(
+        "aatoolbox.datasources.glofas.glofas.cdsapi.Client.retrieve"
+    )
+
+
+@pytest.fixture
+def mock_result(mocker):
+    """Mock the entire Result class."""
+    # This is needed because the reply is changed dynamically
+    # so very difficult to use a Mock object
+    class MockResult:
+        def __init__(self, *args, **kwargs):
+            self.reply = {"state": "completed"}
+            self.state = None
+
+        def update(self):
+            pass
+
+        @mock.create_autospec
+        def download(self, filepath):
+            pass
+
+    return mocker.patch(
+        "aatoolbox.datasources.glofas.glofas.cdsapi.api.Result",
+        return_value=MockResult(),
+    )
+
+
+@pytest.fixture
+def geo_bounding_box():
+    """Input GeoBoundingBox to use."""
+    gbb = GeoBoundingBox(lat_max=1.0, lat_min=-2.2, lon_max=3.3, lon_min=-4.4)
+    return gbb
+
+
+def test_reanalysis_download(
+    mock_country_config,
+    mock_aa_data_dir,
+    mock_retrieve,
+    mock_result,
+    geo_bounding_box,
+):
+    """
+    Test GloFAS reanalysis download.
+
+    Test that the query generated by the download method of GlofasReanlysis
+    with default parameters is as expected
+    """
+    glofas_reanalysis = GlofasReanalysis(
+        country_config=mock_country_config,
+        geo_bounding_box=geo_bounding_box,
+        start_date=date(year=2021, month=1, day=1),
+        end_date=date(year=2021, month=12, day=31),
+    )
+    glofas_reanalysis.download()
+    expected_retrieve_args = {
+        "name": "cems-glofas-historical",
+        "request": {
+            "variable": "river_discharge_in_the_last_24_hours",
+            "format": "grib",
+            "product_type": "consolidated",
+            "system_version": "version_3_1",
+            "hydrological_model": "lisflood",
+            "hyear": "2021",
+            "hmonth": [str(x + 1).zfill(2) for x in range(12)],
+            "hday": [str(x + 1).zfill(2) for x in range(31)],
+            "area": [1.05, -4.45, -2.25, 3.35],
+        },
+    }
+    expected_result_path = Path(
+        f"{mock_aa_data_dir}/public/raw/{mock_country_config.iso3}"
+        f"/glofas/cems-glofas-historical/"
+        f"{mock_country_config.iso3}_"
+        f"cems-glofas-historical_2021_Np1d05Sm2d25Ep3d35Wm4d45.grib"
+    )
+    mock_retrieve.assert_called_with(**expected_retrieve_args)
+    mock_result.return_value.download.assert_called_with(
+        mock.ANY, expected_result_path
+    )
+
+
+def test_forecast_download(
+    mock_country_config,
+    mock_aa_data_dir,
+    mock_retrieve,
+    mock_result,
+    geo_bounding_box,
+):
+    """
+    Test GloFAS forecast download.
+
+    Test that the query generated by the download method of GlofasForecast
+    with default parameters is as expected
+    """
+    glofas_forecast = GlofasForecast(
+        country_config=mock_country_config,
+        geo_bounding_box=geo_bounding_box,
+        leadtime_max=3,
+        start_date=date(year=2022, month=1, day=1),
+        end_date=date(year=2022, month=1, day=1),
+    )
+    glofas_forecast.download()
+    expected_retrieve_args = {
+        "name": "cems-glofas-forecast",
+        "request": {
+            "variable": "river_discharge_in_the_last_24_hours",
+            "format": "grib",
+            "product_type": [
+                "control_forecast",
+                "ensemble_perturbed_forecasts",
+            ],
+            "system_version": "operational",
+            "hydrological_model": "lisflood",
+            "year": "2022",
+            "month": "01",
+            "day": "01",
+            "area": [1.05, -4.45, -2.25, 3.35],
+            "leadtime_hour": ["24", "48", "72"],
+        },
+    }
+
+    expected_result_path = Path(
+        f"{mock_aa_data_dir}/public/raw/{mock_country_config.iso3}/"
+        f"glofas/cems-glofas-forecast/"
+        f"{mock_country_config.iso3}_"
+        f"cems-glofas-forecast_2022-01-01_ltmax03d_Np1d05Sm2d25Ep3d35Wm4d45"
+        f".grib"
+    )
+    mock_retrieve.assert_called_with(**expected_retrieve_args)
+    mock_result.return_value.download.assert_called_with(
+        mock.ANY, expected_result_path
+    )
+
+
+def test_reforecast_download(
+    mock_country_config,
+    mock_aa_data_dir,
+    mock_retrieve,
+    mock_result,
+    geo_bounding_box,
+):
+    """
+    Test GloFAS reforecast download.
+
+    Test that the query generated by the download method of
+    GlofasReforecast with default parameters is as expected
+    """
+    glofas_reforecast = GlofasReforecast(
+        country_config=mock_country_config,
+        geo_bounding_box=geo_bounding_box,
+        leadtime_max=3,
+        start_date=date(year=2018, month=1, day=1),
+        end_date=date(year=2018, month=1, day=31),
+    )
+    glofas_reforecast.download()
+    expected_retrieve_args = {
+        "name": "cems-glofas-reforecast",
+        "request": {
+            "variable": "river_discharge_in_the_last_24_hours",
+            "format": "grib",
+            "product_type": [
+                "control_reforecast",
+                "ensemble_perturbed_reforecasts",
+            ],
+            "system_version": "version_3_1",
+            "hydrological_model": "lisflood",
+            "hyear": "2018",
+            "hmonth": "01",
+            "hday": [str(x + 1).zfill(2) for x in range(31)],
+            "area": [1.05, -4.45, -2.25, 3.35],
+            "leadtime_hour": ["24", "48", "72"],
+        },
+    }
+    expected_result_path = Path(
+        f"{mock_aa_data_dir}/public/raw/{mock_country_config.iso3}/"
+        f"glofas/cems-glofas-reforecast/"
+        f"{mock_country_config.iso3}_"
+        f"cems-glofas-reforecast_2018-01_ltmax03d_Np1d05Sm2d25Ep3d35Wm4d45"
+        f".grib"
+    )
+    mock_retrieve.assert_called_with(**expected_retrieve_args)
+    mock_result.return_value.download.assert_called_with(
+        mock.ANY, expected_result_path
+    )
diff --git a/tests/datasources/test_glofas_process.py b/tests/datasources/test_glofas_process.py
new file mode 100644
index 00000000..292f3581
--- /dev/null
+++ b/tests/datasources/test_glofas_process.py
@@ -0,0 +1,274 @@
+"""Tests for GloFAS processing."""
+from datetime import date
+from typing import List, Union
+
+import numpy as np
+import pandas as pd
+import pytest
+import xarray as xr
+
+from aatoolbox import (
+    CountryConfig,
+    GeoBoundingBox,
+    GlofasForecast,
+    GlofasReanalysis,
+    GlofasReforecast,
+)
+
+
+class TestProcess:
+    """Tests for GloFAS processing."""
+
+    geo_bounding_box = GeoBoundingBox(
+        lat_max=1, lat_min=-2, lon_max=3, lon_min=-4
+    )
+    numbers = [0, 1, 2, 3, 4, 5, 6]
+
+    def get_raw_data(
+        self,
+        number_coord: Union[List[int], int] = None,
+        include_step: bool = False,
+        include_history: bool = False,
+        dis24: np.ndarray = None,
+        single_day=False,
+    ) -> xr.Dataset:
+        """
+        Construct a simple fake GloFAS xarray dataset.
+
+        Parameters
+        ----------
+        number_coord : list or int, default = None
+            The ensemble number coordinate
+        include_step :  bool, default = False
+            Whether to include the forecast step coordinate
+        include_history : bool, default = False
+            Whether to include the history attribute
+        dis24 : np.ndarray, default = None
+            Optional array of discharge values, that should have the combined
+            dimensions of the coordinates. If not passed, generated using
+            random numbers.
+        single_day: boolean, default = False
+            Only create a single day of data, required for the forecasts
+
+        Returns
+        -------
+        Simplified GloFAS xarray dataset
+        """
+        rng = np.random.default_rng(12345)
+        coords = {}
+        if number_coord is not None:
+            coords["number"] = number_coord
+        if single_day:
+            coords["time"] = np.datetime64("2014-09-06")
+        else:
+            coords["time"] = pd.date_range("2014-09-06", periods=2)
+        if include_step:
+            coords["step"] = [np.datetime64(n + 1, "D") for n in range(5)]
+        coords["latitude"] = (
+            np.arange(
+                start=self.geo_bounding_box.lat_min,
+                stop=self.geo_bounding_box.lat_max + 2,
+                step=0.1,
+            )
+            - 0.05
+        )
+        coords["longitude"] = (
+            np.arange(
+                start=self.geo_bounding_box.lon_min,
+                stop=self.geo_bounding_box.lon_max + 2,
+                step=0.1,
+            )
+            - 0.05
+        )
+        dims = list(coords.keys())
+        if number_coord is not None and isinstance(number_coord, int):
+            dims.remove("number")
+        if single_day:
+            dims.remove("time")
+        if dis24 is None:
+            dis24 = 5000 + 100 * rng.random([len(coords[dim]) for dim in dims])
+        attrs = {}
+        if include_history:
+            attrs = {"history": "fake history"}
+        return xr.Dataset({"dis24": (dims, dis24)}, coords=coords, attrs=attrs)
+
+    @pytest.fixture()
+    def mock_ensemble_raw(self):
+        """
+        Create fake raw ensemble data.
+
+        For the forecast and reforecast, generate the raw data, which consists
+        of the control and perturbed forecast, and combine the discharge
+        values for the processed data. Return both xarray datasets and the
+        combined array.
+        """
+
+        def _mock_ensemble_raw(
+            single_day: bool = False,
+        ) -> (xr.Dataset, xr.Dataset, np.ndarray):
+            cf_raw = self.get_raw_data(
+                number_coord=self.numbers[0],
+                include_step=True,
+                include_history=True,
+                single_day=single_day,
+            )
+            pf_raw = self.get_raw_data(
+                number_coord=self.numbers[1:],
+                include_step=True,
+                include_history=True,
+                single_day=single_day,
+            )
+            expected_dis24 = np.concatenate(
+                (
+                    cf_raw["dis24"].values[np.newaxis, ...],
+                    pf_raw["dis24"].values,
+                )
+            )
+            return cf_raw, pf_raw, expected_dis24
+
+        return _mock_ensemble_raw
+
+    def get_processed_data(
+        self,
+        country_config: CountryConfig,
+        number_coord: [List[int], int] = None,
+        include_step: bool = False,
+        dis24: np.ndarray = None,
+        single_day: bool = False,
+    ) -> xr.Dataset:
+        """
+        Create a simplified fake processed GloFAS dataset.
+
+        Parameters
+        ----------
+        country_config : CountryConfig
+            Country configuration object
+        number_coord : list or int, default = None
+            The ensemble number coordinate
+        include_step : bool, default = False
+            Whether to include the forecast step coordinate
+        dis24 : np.ndarray, default = None
+            Optional array of discharge values, that should have the combined
+            dimensions of the coordinates. If not passed, generated using
+            random numbers.
+        single_day: boolean, default = False
+            Only create a single day of data, required for the forecasts
+
+        Returns
+        -------
+        GloFAS processed xarray dataset
+        """
+        raw_data = self.get_raw_data(
+            number_coord=number_coord,
+            include_step=include_step,
+            dis24=dis24,
+            single_day=single_day,
+        )
+        coords = {}
+        if number_coord is not None:
+            coords = {"number": number_coord}
+        if not single_day:
+            coords["time"] = raw_data.time
+        if include_step:
+            coords["step"] = raw_data.step
+        return xr.Dataset(
+            {
+                reporting_point.name: (
+                    list(coords.keys()),
+                    raw_data["dis24"]
+                    .sel(
+                        longitude=reporting_point.lon,
+                        latitude=reporting_point.lat,
+                        method="nearest",
+                    )
+                    .data,
+                )
+                for reporting_point in country_config.glofas.reporting_points
+            },
+            coords=coords,
+        )
+
+    @pytest.fixture()
+    def mock_processed_data_reanalysis(
+        self, mocker, mock_country_config
+    ) -> xr.Dataset:
+        """Create fake processed GloFAS reanalysis data."""
+        mocker.patch(
+            "aatoolbox.datasources.glofas.reanalysis.xr.load_dataset",
+            return_value=self.get_raw_data(),
+        )
+        return self.get_processed_data(country_config=mock_country_config)
+
+    @pytest.fixture()
+    def mock_processed_data_forecast(
+        self, mock_ensemble_raw, mocker, mock_country_config
+    ):
+        """Create fake processed GloFAS forecast or reforecast data."""
+
+        def _mock_processed_dta_forecast(single_day: bool = False):
+            cf_raw, pf_raw, expected_dis24 = mock_ensemble_raw(
+                single_day=single_day
+            )
+            mocker.patch(
+                "aatoolbox.datasources.glofas.forecast.xr.load_dataset",
+                side_effect=[cf_raw, pf_raw],
+            )
+            return self.get_processed_data(
+                country_config=mock_country_config,
+                number_coord=self.numbers,
+                include_step=True,
+                dis24=mock_ensemble_raw(single_day=single_day)[2],
+                single_day=single_day,
+            )
+
+        return _mock_processed_dta_forecast
+
+    def test_reanalysis_process(
+        self, mock_country_config, mock_processed_data_reanalysis
+    ):
+        """Test GloFAS reanalysis process method."""
+        glofas_reanalysis = GlofasReanalysis(
+            country_config=mock_country_config,
+            geo_bounding_box=self.geo_bounding_box,
+            start_date=date(year=2021, month=1, day=1),
+            end_date=date(year=2021, month=12, day=31),
+        )
+        output_filepath = glofas_reanalysis.process()[0]
+        # use open_dataset since load_dataset is patched
+        with xr.open_dataset(output_filepath) as output_ds:
+            assert output_ds.equals(mock_processed_data_reanalysis)
+
+    def test_reforecast_process(
+        self, mock_country_config, mock_processed_data_forecast
+    ):
+        """Test GloFAS reforecast process method."""
+        target_dataset = mock_processed_data_forecast()
+        glofas_reforecast = GlofasReforecast(
+            country_config=mock_country_config,
+            geo_bounding_box=self.geo_bounding_box,
+            leadtime_max=3,
+            start_date=date(year=2018, month=1, day=1),
+            end_date=date(year=2018, month=1, day=31),
+        )
+        output_filepath = glofas_reforecast.process()[0]
+        with xr.open_dataset(output_filepath) as output_ds:
+            assert output_ds.equals(target_dataset)
+
+    def test_forecast_process(
+        self,
+        mock_country_config,
+        mock_processed_data_forecast,
+        mock_ensemble_raw,
+    ):
+        """Test GloFAS forecast process method."""
+        target_dataset = mock_processed_data_forecast(single_day=True)
+        glofas_forecast = GlofasForecast(
+            country_config=mock_country_config,
+            geo_bounding_box=self.geo_bounding_box,
+            leadtime_max=3,
+            start_date=date(year=2022, month=1, day=1),
+            end_date=date(year=2022, month=1, day=1),
+        )
+        output_filepath = glofas_forecast.process()[0]
+        with xr.open_dataset(output_filepath) as output_ds:
+            assert output_ds.equals(target_dataset)
diff --git a/tests/datasources/test_iri.py b/tests/datasources/test_iri.py
index 9ba1d574..22517412 100644
--- a/tests/datasources/test_iri.py
+++ b/tests/datasources/test_iri.py
@@ -1,10 +1,8 @@
 """Tests for the IRI module."""
-
-from pathlib import Path
-
 import cftime
 import numpy as np
 import pytest
+import requests
 import xarray as xr
 from xarray.coding.cftimeindex import CFTimeIndex
 
@@ -108,6 +106,29 @@ def test_download_call_dominant(
     )
 
 
+@pytest.fixture
+def mock_requests(mocker):
+    """Mock requests in the download function."""
+    requests_mock = mocker.patch(
+        "aatoolbox.datasources.iri.iri_seasonal_forecast.requests.get"
+    )
+    mocker.patch.dict(
+        "aatoolbox.datasources.iri.iri_seasonal_forecast.os.environ",
+        {"IRI_AUTH": FAKE_IRI_AUTH},
+    )
+
+    return requests_mock
+
+
+def test_download_wrong_auth(
+    mock_iri, mock_requests, mock_aa_data_dir, mock_country_config
+):
+    """Check that wrong download headers raise an error."""
+    iri = mock_iri(prob_forecast=True)
+    with pytest.raises(requests.RequestException):
+        iri.download()
+
+
 def test_process(mocker, mock_iri, mock_aa_data_dir, mock_country_config):
     """Test process for IRI forecast."""
     ds = xr.DataArray(
@@ -161,7 +182,7 @@ def test_process_if_download_not_called(mock_iri):
     iri = mock_iri()
     # Make sure file doesn't exist
     if iri._get_raw_path().exists():
-        Path.unlink(iri._get_raw_path())
+        iri._get_raw_path().unlink()
     with pytest.raises(FileNotFoundError) as excinfo:
         iri.process()
     assert (
@@ -213,7 +234,7 @@ def test_load_if_process_not_called(mock_iri):
     iri = mock_iri()
     # Make sure file doesn't exist
     if iri._get_processed_path().exists():
-        Path.unlink(iri._get_processed_path())
+        iri._get_processed_path().unlink()
     with pytest.raises(FileNotFoundError) as excinfo:
         iri.load()
     assert (
diff --git a/tests/datasources/test_ndvi.py b/tests/datasources/test_ndvi.py
new file mode 100644
index 00000000..f3043956
--- /dev/null
+++ b/tests/datasources/test_ndvi.py
@@ -0,0 +1,303 @@
+"""Tests for the USGS NDVI module."""
+import tempfile
+from datetime import date, datetime, timedelta
+from pathlib import Path
+
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import pytest
+import xarray as xr
+from shapely.geometry import Polygon
+
+from aatoolbox import (
+    UsgsNdviMedianAnomaly,
+    UsgsNdviPctMedian,
+    UsgsNdviSmoothed,
+    UsgsNdviYearDifference,
+)
+
+DATASOURCE_BASE_DIR = "usgs_ndvi"
+
+
+@pytest.fixture
+def mock_ndvi(mock_country_config):
+    """Create USGS NDVI class with mock country config."""
+    start_date = (2019, 36)
+    end_date = (2020, 2)
+
+    instantiator = {
+        "smoothed": UsgsNdviSmoothed,
+        "pct_median": UsgsNdviPctMedian,
+        "anomaly": UsgsNdviMedianAnomaly,
+        "difference": UsgsNdviYearDifference,
+    }
+
+    def _mock_ndvi(variable: str = "smoothed"):
+
+        ndvi = instantiator[variable](
+            country_config=mock_country_config,
+            start_date=start_date,
+            end_date=end_date,
+        )
+        return ndvi
+
+    return _mock_ndvi
+
+
+@pytest.fixture
+def mock_download(mocker, mock_ndvi):
+    """
+    Call download with mocked download.
+
+    Return number of calls to the internal
+    download method and the returned filepath.
+    """
+    download_mock = mocker.patch(
+        (
+            "aatoolbox.datasources.usgs."
+            "ndvi_base._UsgsNdvi._download_ndvi_dekad"
+        )
+    )
+
+    def _mock_download(variable: str = "smoothed"):
+        ndvi = mock_ndvi(variable=variable)
+        fp = ndvi.download()
+        download_calls = download_mock.call_count
+        return fp, download_calls
+
+    return _mock_download
+
+
+@pytest.fixture
+def mock_determine_process_dates(mocker, mock_ndvi):
+    """
+    Call _determine_process_dates() with mocked _load().
+
+    Return number of calls to the internal
+    download method and the returned filepath.
+    """
+    mocker.patch(
+        ("aatoolbox.datasources.usgs.ndvi_base._UsgsNdvi._load"),
+        return_value=pd.DataFrame(
+            {
+                "year": [2019, 2020, 2020],
+                "dekad": [36, 1, 2],
+                "modified": [datetime.now() + timedelta(days=10)] * 3,
+            }
+        ),
+    )
+
+    mocker.patch(
+        "aatoolbox.datasources.usgs.ndvi_base._UsgsNdvi._get_raw_path",
+        return_value=Path(tempfile.mkstemp()[1]),
+    )
+
+    def _mock_determine_process_dates(clobber: bool, dates_to_process: list):
+        ndvi = mock_ndvi(variable="smoothed")
+        return ndvi._determine_process_dates(
+            clobber=clobber,
+            filepath="fake_filepath",
+            dates_to_process=dates_to_process,
+        )
+
+    return _mock_determine_process_dates
+
+
+@pytest.fixture
+def gdf():
+    """Create geodataframe for testing."""
+    d = {
+        "name": ["area_a", "area_b"],
+        "geometry": [
+            Polygon([(0, 0), (0, 2), (2, 2), (2, 0)]),
+            Polygon([(2, 0), (2, 2), (3, 2), (3, 0)]),
+        ],
+    }
+    return gpd.GeoDataFrame(d)
+
+
+@pytest.fixture
+def da():
+    """
+    Create raster input with date coordinate.
+
+    Used as side effect later to create 3
+    consecutive raster values for processing.
+    """
+
+    def _da(year_dekad):
+        da = xr.DataArray(
+            [[[1 + year_dekad[1], 2, 3], [4, 5 + year_dekad[1], 6]]],
+            dims=("band", "y", "x"),
+            coords={
+                "y": [1.5, 0.5],
+                "x": [0.5, 1.5, 2.5],
+            },
+        ).rio.write_crs("EPSG:4326")
+        return da
+
+    return _da
+
+
+def test_download_smoothed(mock_download, mock_aa_data_dir):
+    """Test download for all NDVI variables."""
+    fp, call_count = mock_download(variable="smoothed")
+
+    assert fp == (mock_aa_data_dir / f"public/raw/glb/{DATASOURCE_BASE_DIR}")
+    assert call_count == 3
+
+
+def test_download_pct_median(mock_download, mock_aa_data_dir):
+    """Test download for all NDVI variables."""
+    fp, call_count = mock_download(variable="pct_median")
+
+    assert fp == (mock_aa_data_dir / f"public/raw/glb/{DATASOURCE_BASE_DIR}")
+
+    assert call_count
+
+
+def test_download_anomaly(mock_download, mock_aa_data_dir):
+    """Test download for all NDVI variables."""
+    fp, call_count = mock_download(variable="anomaly")
+
+    assert fp == (mock_aa_data_dir / f"public/raw/glb/{DATASOURCE_BASE_DIR}")
+
+    assert call_count
+
+
+def test_download_difference(mock_download, mock_aa_data_dir):
+    """Test download for all NDVI variables."""
+    fp, call_count = mock_download(variable="difference")
+
+    assert fp == (mock_aa_data_dir / f"public/raw/glb/{DATASOURCE_BASE_DIR}")
+
+    assert call_count
+
+
+def test_process_and_load(
+    mocker, mock_ndvi, mock_aa_data_dir, mock_country_config, gdf, da
+):
+    """Test processing NDVI values."""
+    ndvi = mock_ndvi(variable="smoothed")
+
+    # TODO: now created `load_raw` to be able to mock but would like
+    #  to do it from xr.load_dataset directly
+    mocker.patch(
+        "aatoolbox.datasources.usgs.ndvi_base.rioxarray.open_rasterio",
+        side_effect=[da((2019, 36)), da((2020, 1)), da((2020, 2))] * 6,
+    )
+
+    mocker.patch(
+        "aatoolbox.datasources.usgs.ndvi_base._UsgsNdvi._get_raw_path",
+        return_value=Path(tempfile.mkstemp()[1]),
+    )
+
+    processed_dir = ndvi.process(gdf=gdf, feature_col="name")
+    assert processed_dir == (
+        mock_aa_data_dir / f"public/processed/{mock_country_config.iso3}/"
+        f"{DATASOURCE_BASE_DIR}"
+    )
+
+    df_processed = ndvi.load(feature_col="name")
+
+    expected_dates = [
+        date(year=2019, month=12, day=21),
+        date(year=2020, month=1, day=1),
+        date(year=2020, month=1, day=11),
+    ]
+
+    expected_dates = [d for ds in expected_dates for d in [ds] * 2]
+
+    assert np.array_equal(df_processed.name.values, ["area_a", "area_b"] * 3)
+    assert np.array_equal(
+        pd.to_datetime(df_processed.date).dt.date, expected_dates
+    )
+    assert np.array_equal(df_processed.dekad.values, [36, 36, 1, 1, 2, 2])
+    assert np.array_equal(df_processed.loc[:, "count"].values, [4, 2] * 3)
+    assert np.array_equal(
+        df_processed.loc[:, "min"].values, [2, 3, 2, 3, 2, 3]
+    )
+    assert np.allclose(
+        df_processed.loc[:, "std"].values,
+        [18.069311, 1.5, 1.6583124, 1.5, 1.87082869, 1.5],
+    )
+    assert np.array_equal(
+        df_processed.loc[:, "mean"].values, [21.0, 4.5, 3.5, 4.5, 4.0, 4.5]
+    )
+    assert np.array_equal(
+        df_processed.loc[:, "max"].values, [41, 6, 6, 6, 7, 6]
+    )
+
+
+def test_process_if_download_not_called(mock_ndvi, gdf):
+    """Test that correct error message raised."""
+    ndvi = mock_ndvi(variable="smoothed")
+    with pytest.raises(FileNotFoundError):
+        ndvi.process(gdf=gdf, feature_col="name")
+
+
+def test_load_if_process_not_called(mock_ndvi):
+    """Test that correct error message raised."""
+    ndvi = mock_ndvi(variable="smoothed")
+    with pytest.raises(FileNotFoundError):
+        ndvi.load(feature_col="name")
+
+
+def test_fp_year_dekad():
+    """Test that year and dekad extracted from FP."""
+    fp = Path("ea2022_10pct.tif")
+    assert UsgsNdviPctMedian._fp_year_dekad(fp) == [2022, 10]
+
+
+def test_get_url(mock_ndvi):
+    """Test that URL generated correctly."""
+    fn = "fake_filename"
+    ndvi = mock_ndvi(variable="pct_median")
+    expected_fn = (
+        f"https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/fews/web/"
+        f"{ndvi._datasource_config.area_url}/dekadal/emodis"
+        f"/ndvi_c6/{ndvi._data_variable_url}/"
+        f"downloads/dekadal/{fn}.zip"
+    )
+    assert ndvi._get_url(fn) == expected_fn
+
+
+def test_process_dates_clobber_true(mock_determine_process_dates):
+    """Test process dates clobbers properly."""
+    dps = [(2019, 36), (2020, 1), (2020, 2)]
+    test_dps, df = mock_determine_process_dates(
+        clobber=True, dates_to_process=dps
+    )
+    assert len(test_dps) == 3
+    assert df.size == 0
+
+
+def test_process_dates_clobber_false(mock_determine_process_dates):
+    """Test process dates doesn't clobber unnecessarily."""
+    dps = [(2019, 36), (2020, 1), (2020, 2)]
+    test_dps, df = mock_determine_process_dates(
+        clobber=False, dates_to_process=dps
+    )
+    assert len(test_dps) == 0
+    assert df.shape[0] == 3
+
+
+def test_process_dates_clobber_true_additional(mock_determine_process_dates):
+    """Test process dates clobbers properly."""
+    dps = [(2019, 35), (2019, 36), (2020, 1), (2020, 2), (2020, 3)]
+    test_dps, df = mock_determine_process_dates(
+        clobber=True, dates_to_process=dps
+    )
+    assert len(test_dps) == 5
+    assert df.size == 0
+
+
+def test_process_dates_clobber_false_additional(mock_determine_process_dates):
+    """Test process dates doesn't clobber unnecessarily."""
+    dps = [(2019, 35), (2019, 36), (2020, 1), (2020, 2), (2020, 3)]
+    test_dps, df = mock_determine_process_dates(
+        clobber=False, dates_to_process=dps
+    )
+    assert len(test_dps) == 2
+    assert df.shape[0] == 3
diff --git a/tests/test_country_config.py b/tests/test_country_config.py
index 0174491d..4a54c309 100644
--- a/tests/test_country_config.py
+++ b/tests/test_country_config.py
@@ -72,7 +72,7 @@ def test_uppercase_config_iso3(mock_parse_yaml):
     assert country_config.iso3 == "abc"
 
 
-def test_validate_codab_layer_base_name(mock_parse_yaml):
+def test_codab_validate_layer_base_name(mock_parse_yaml):
     """Test that layer basename requires correct placeholder."""
     config_base = {
         "iso3": "abc",
@@ -94,7 +94,35 @@ def test_validate_codab_layer_base_name(mock_parse_yaml):
         create_country_config(iso3="abc")
 
 
-def test_validate_fewsnet_region_name(mock_parse_yaml):
+def test_codab_validate_admin_level_max(mock_parse_yaml):
+    """Test that admin level can only be between 0 and 4."""
+    config_base = {
+        "iso3": "abc",
+        "codab": {
+            "hdx_dataset_name": "fake_dataset_name",
+            "layer_base_name": "layer_base_name_{admin_level}",
+            "admin_level_max": 1,
+        },
+    }
+    config_correct = copy.deepcopy(config_base)
+    config_incorrect_a = copy.deepcopy(config_base)
+    config_incorrect_a["codab"]["admin_level_max"] = -1
+    config_incorrect_b = copy.deepcopy(config_base)
+    config_incorrect_b["codab"]["admin_level_max"] = 5
+    mock_parse_yaml(
+        output_list=[config_correct, config_incorrect_a, config_incorrect_b]
+    )
+    # Check that correct config runs without issue
+    create_country_config(iso3="abc")
+    # Check that incorrect config raises error for config incorrect a
+    with pytest.raises(ValueError):
+        create_country_config(iso3="abc")
+    # Check that incorrect config raises error for config incorrect b
+    with pytest.raises(ValueError):
+        create_country_config(iso3="abc")
+
+
+def test_fewsnet_validate_region_name(mock_parse_yaml):
     """Test that fewsnet requires correct region name."""
     config_base = {
         "iso3": "abc",
diff --git a/tests/test_check_file_existence.py b/tests/utils/test_check_file_existence.py
similarity index 55%
rename from tests/test_check_file_existence.py
rename to tests/utils/test_check_file_existence.py
index 94763f19..255c8015 100644
--- a/tests/test_check_file_existence.py
+++ b/tests/utils/test_check_file_existence.py
@@ -1,5 +1,7 @@
 """Tests for check_file_existence decorator."""
 
+import logging
+
 import pytest
 
 from aatoolbox.utils.check_file_existence import check_file_existence
@@ -11,24 +13,38 @@ def downloader(filepath, clobber):
     return "a"
 
 
-def test_fp_exists_no_clobber(tmp_path):
+def test_fp_exists_no_clobber(tmp_path, caplog):
     """Test that filepath returned if Path exists and clobber False."""
-    tmp_path.touch()
+    caplog.set_level(logging.INFO)
     output_filepath = downloader(filepath=tmp_path, clobber=False)
     assert output_filepath == tmp_path
+    assert (
+        f"File {tmp_path} exists and clobber set to "
+        "False, using existing file." in caplog.text
+    )
 
 
-def test_fp_exists_clobber(tmp_path):
+def test_fp_exists_clobber(tmp_path, caplog):
     """Test that function returned if Path exists and clobber True."""
-    tmp_path.touch()
+    caplog.set_level(logging.INFO)
     output_filepath = downloader(filepath=tmp_path, clobber=True)
     assert output_filepath == "a"
+    assert (
+        f"File {tmp_path} exists and clobber set to "
+        "True, overwriting existing file." in caplog.text
+    )
 
 
-def test_fp_not_exists(tmp_path):
+def test_fp_not_exists(tmp_path, caplog):
     """Test that function returned if Path does not exist."""
-    output_filepath = downloader(filepath=tmp_path, clobber=True)
+    caplog.set_level(logging.INFO)
+    path_new = tmp_path / "new_path"
+    output_filepath = downloader(filepath=path_new, clobber=False)
     assert output_filepath == "a"
+    assert (
+        f"File {path_new} does not exist and clobber set to "
+        "False, downloading new file." in caplog.text
+    )
 
 
 def test_key_error(tmp_path):
diff --git a/tests/utils/test_dates.py b/tests/utils/test_dates.py
new file mode 100644
index 00000000..e859be92
--- /dev/null
+++ b/tests/utils/test_dates.py
@@ -0,0 +1,89 @@
+"""Test dates module."""
+
+from datetime import date, datetime
+
+import pytest
+
+from aatoolbox.utils import dates
+
+
+def test_get_date():
+    """Tests getting date from string."""
+    desired_date = date(year=2018, month=3, day=14)
+    desired_str = "2018-03-14"
+    desired_datetime = datetime(year=2018, month=3, day=14)
+    assert dates.get_date_from_user_input(desired_str) == desired_date
+    assert dates.get_date_from_user_input(desired_date) == desired_date
+    assert dates.get_date_from_user_input(desired_datetime) == desired_datetime
+
+
+def test_get_date_value_error():
+    """Tests getting date errors."""
+    with pytest.raises(ValueError):
+        dates.get_date_from_user_input("2013-3-14")
+    with pytest.raises(TypeError):
+        dates.get_date_from_user_input(None)
+
+
+def test_get_dekadal_date():
+    """Tests get dekadal date."""
+    desired_dekad = (2013, 8)
+    desired_str = "2013-03-14"
+    assert dates.get_dekadal_date(desired_str) == desired_dekad
+    assert dates.get_dekadal_date(desired_dekad) == desired_dekad
+    assert (
+        dates.get_dekadal_date(date.fromisoformat(desired_str))
+        == desired_dekad
+    )
+    assert (
+        dates.get_dekadal_date(input_date=None, default_date=desired_dekad)
+        == desired_dekad
+    )
+
+
+def test_get_dekadal_date_value_error():
+    """Test that value error raised for bad inputs."""
+    with pytest.raises(ValueError):
+        dates.get_dekadal_date("2020-01-1")
+    with pytest.raises(ValueError):
+        dates.get_dekadal_date((20, 17))
+
+
+def test_dekad_to_date():
+    """Test conversion from dekad to date."""
+    assert dates.dekad_to_date((2020, 1)) == date.fromisoformat("2020-01-01")
+    assert dates.dekad_to_date((2016, 20)) == date.fromisoformat("2016-07-11")
+    assert dates.dekad_to_date((2024, 36)) == date.fromisoformat("2024-12-21")
+
+
+def test_date_to_dekad():
+    """Test conversion from dekad to date."""
+    assert dates.date_to_dekad(date.fromisoformat("2020-01-09")) == (2020, 1)
+    assert dates.date_to_dekad(date.fromisoformat("2016-07-20")) == (2016, 20)
+    assert dates.date_to_dekad(date.fromisoformat("2024-12-31")) == (2024, 36)
+
+
+def test_compare_dekads():
+    """Test comparing dekads."""
+    dekad1 = (2019, 36)
+    dekad2 = (2020, 1)
+    assert dates.compare_dekads_gt(dekad2, dekad1)
+    assert dates.compare_dekads_gte(dekad2, dekad1)
+    assert dates.compare_dekads_lt(dekad1, dekad2)
+    assert dates.compare_dekads_lte(dekad1, dekad2)
+    assert dates.compare_dekads_gte(dekad2, dekad2)
+    assert dates.compare_dekads_lte(dekad2, dekad2)
+
+
+def test_expand_dekads():
+    """Test expanding dekads."""
+    dekad1 = (2019, 33)
+    dekad2 = (2020, 3)
+    dekads = dates.expand_dekads(dekad1, dekad2)
+    assert len(dekads) == 7
+    assert dekads[0] == dekad1
+    assert dekads[6] == dekad2
+    assert dekads[4] == (2020, 1)
+    # error raised for reverse attempt
+    with pytest.raises(ValueError):
+        dates.expand_dekads(dekad2, dekad1)
diff --git a/tests/test_geoboundingbox.py b/tests/utils/test_geoboundingbox.py
similarity index 100%
rename from tests/test_geoboundingbox.py
rename to tests/utils/test_geoboundingbox.py
diff --git a/tests/test_hdx_api.py b/tests/utils/test_hdx_api.py
similarity index 100%
rename from tests/test_hdx_api.py
rename to tests/utils/test_hdx_api.py
diff --git a/tests/test_raster.py b/tests/utils/test_raster.py
similarity index 83%
rename from tests/test_raster.py
rename to tests/utils/test_raster.py
index db4b45c4..daf395d3 100644
--- a/tests/test_raster.py
+++ b/tests/utils/test_raster.py
@@ -12,11 +12,28 @@
 from rioxarray.exceptions import DimensionError, MissingCRS
 from shapely.geometry import Polygon
 
-from aatoolbox.utils import raster
+import aatoolbox
 
-suite = unittest.TestSuite()
-suite.addTest(doctest.DocTestSuite(raster))
-unittest.TextTestRunner().run(suite)
+
+def test_doctest_suite():
+    """
+    Test docstrings in raster module.
+
+    Checks if there are any failures in the doctest
+    suite, meaning that there is some error within
+    the docstrings of the raster module. This is used
+    to avoid running doctest on all docstrings in the
+    library, which would require many exceptions given
+    the extent of downloading and processing other modules
+    rely on.
+
+    Method found in:
+    https://vladyslav-krylasov.medium.com/discover-unit-tests-and-doctests-in-one-run-c5504aea86bd
+    """
+    suite = unittest.TestSuite()
+    suite.addTest(doctest.DocTestSuite(aatoolbox.utils.raster))
+    runner = unittest.TextTestRunner(verbosity=2).run(suite)
+    assert not runner.failures
 
 
 @pytest.fixture
@@ -83,12 +100,12 @@ def expected_2d():
     """Create expected comped stats for 2d raster."""
     df = pd.DataFrame(
         {
-            "mean_name": {0: 3.0, 1: 4.5},
-            "std_name": {0: 1.5811388300841898, 1: 1.5},
-            "min_name": {0: 1, 1: 3},
-            "max_name": {0: 5, 1: 6},
-            "sum_name": {0: 12.0, 1: 9.0},
-            "count_name": {0: 4, 1: 2},
+            "mean": {0: 3.0, 1: 4.5},
+            "std": {0: 1.5811388300841898, 1: 1.5},
+            "min": {0: 1, 1: 3},
+            "max": {0: 5, 1: 6},
+            "sum": {0: 12.0, 1: 9.0},
+            "count": {0: 4, 1: 2},
             "name": {0: "area_a", 1: "area_b"},
         },
     )
@@ -108,7 +125,7 @@ def test_compute_raster_stats_2d(da_2d, gdf, expected_2d):
     """Compute raster stats without time dimensions."""
     result = da_2d.aat.compute_raster_stats(gdf, "name")
     assert_frame_equal(result, expected_2d, check_dtype=False)
-    expected_2d.insert(loc=6, column="95quant_name", value=[4.85, 5.85])
+    expected_2d.insert(loc=6, column="95quant", value=[4.85, 5.85])
     result_pct = da_2d.aat.compute_raster_stats(
         gdf=gdf, feature_col="name", percentile_list=[95]
     )
@@ -134,7 +151,7 @@ def test_compute_raster_stats_3d(ds_3d, gdf, expected_3d):
 def test_compute_raster_stats_da_assertions(da_3d, gdf):
     """Ensure error assertions working in compute raster stats."""
     with pytest.raises(MissingCRS):
-        da_3d.aat._crs = False
+        da_3d.rio._crs = False
         da_3d.aat.compute_raster_stats(gdf=gdf, feature_col="name")