Repository: OuyangWenyu/hydrodataset (language: python)

Issue:
Dev
DAYMET, ERA5-LAND, PMLV2, MOD16A2V105, NLDAS


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index a2609d9..a725cfe 100644
--- a/README.md
+++ b/README.md
@@ -4,6 +4,7 @@ Data downloader and processor for Hydrologic Modeling
 
 ## Data source zoo list
 
+- CAMELS/MOPEX/LAMAH
 - Daymet
 - ECMWF
 - MODIS
@@ -11,27 +12,98 @@ Data downloader and processor for Hydrologic Modeling
 
 More details are shown in the following sections.
 
-## Daymet
+## CAMELS/MOPEX/LAMAH
+
+The CAMELS series data include:
+
+- CAMELS-AUS
+  ([CAMELS-AUS: Hydrometeorological time series and landscape attributes for 222 catchments in Australia](https://essd.copernicus.org/preprints/essd-2020-228/))
+- CAMELS-BR
+  ([CAMELS-BR: Hydrometeorological time series and landscape attributes for 897 catchments in Brazil - link to files](https://doi.org/10.5194/essd-12-2075-2020))
+- CAMELS-CL
+  ([The CAMELS-CL dataset: catchment attributes and meteorology for large sample studies – Chile dataset](https://doi.org/10.5194/hess-22-5817-2018))
+- CAMELS-GB
+  ([CAMELS-GB: Hydrometeorological time series and landscape attributes for 671 catchments in Great Britain](https://doi.org/10.5194/essd-2020-49))
+- CAMELS-US
+  ([The CAMELS data set: catchment attributes and meteorology for large-sample studies](https://doi.org/10.5194/hess-21-5293-2017))
+- CAMELS-YR
+  ([Catchment attributes and meteorology for large sample study in contiguous China](https://doi.org/10.5194/essd-2021-71))
 
-We download and process Daymet data for the 671 basins in [CAMELS](https://ral.ucar.edu/solutions/products/camels).
+We also support [CANOPEX](https://doi.org/10.1002/hyp.10880) (Canada's MOPEX dataset)
+and [LamaH-CE](https://doi.org/10.5194/essd-13-4529-2021) (similar with CAMELS and it is for Central Europe), because we
+use these datasets just like CAMELS, we write similar code in data/data_camels.py.
 
 If you can read Chinese, [this blog](https://github.com/OuyangWenyu/aqualord/blob/master/CAMELS/CAMELS.md) may be a
-quick start for CAMELS.
+quick start for CAMELS (CAMELS-US)
+and [this](https://github.com/OuyangWenyu/aqualord/blob/master/CAMELS/CAMELS-other.md)
+for other CAMELS datasets.
 
-### Downloading the CAMELS dataset
+### Download CAMELS/MOPEX/LAMAH datasets
 
-You can download CAMELS manually from https://ral.ucar.edu/solutions/products/camels ; or you can use the following
-code:
+We recommend downloading the datasets manually, the downloading address are as follows:
+
+- [Download CAMELS-AUS](https://doi.pangaea.de/10.1594/PANGAEA.921850)
+- [Download CAMELS-BR](https://doi.org/10.5281/zenodo.3709337)
+- [Download CAMELS-CL](https://doi.pangaea.de/10.1594/PANGAEA.894885)
+- [Download CAMELS-GB](https://doi.org/10.5285/8344e4f3-d2ea-44f5-8afa-86d2987543a9)
+- [Download CAMELS-US](https://ral.ucar.edu/solutions/products/camels)
+- [Download CAMELS-YR](http://doi.org/10.5281/zenodo.4704017)
+- [Download CANOPEX](http://canopex.etsmtl.net/)
+- [Download LamaH-CE](https://zenodo.org/record/5153305#.YYdEgGBByUk)
+
+You can also use the following code to download CAMELS-US (notice: the unzipped file is 10+ GB):
 
 ```Python
 import os
 import definitions
 from hydrobench.data.data_camels import Camels
 
-camels_path = os.path.join(definitions.DATASET_DIR, "camels")
+# DATASET_DIR is defined in the definitions.py file
+camels_path = os.path.join(definitions.DATASET_DIR, "camels", "camels_us")
 camels = Camels(camels_path, download=True)
 ```
 
+For CAMELS_YR, it is enough to
+download [9_Normal_Camels_YR.zip](https://zenodo.org/record/4704017/files/9_Normal_Camels_YR.zip?download=1)
+
+To download CANOPEX, you have to deal with the GFW. In addtion, there is no attributes data in CANOPEX, we choose an
+alternative: [attributes data](https://osf.io/7fn4c/) from [HYSETS](https://doi.org/10.1038/s41597-020-00583-2)
+
+After downloading, puteach dataset in one directory, the following file-organization is recommended:
+
+```Directory
+camels
+│
+└── camels_aus
+    └── 01_id_name_metadata.zip
+    └── 02_location_boundary_area.zip
+    └── ...
+└── camels_br
+    └── ...
+└── camels_cl
+    └── ...
+└── camels_gb
+    └── ... 
+└── camels_us
+    └── ... 
+└── camels_yr
+    └── ... 
+canopex
+    └── Boundaries.zip
+    └── ...   
+lamah_ce
+    └── 2_LamaH-CE_daily.tar.gz
+    └── ...   
+```
+
+### Process datasets
+
+All methods for processing CAMELS datasets are written in Camels class in hydrobench/data/data_camels.py.
+
+## Daymet
+
+We download and process Daymet data for the 671 basins in [CAMELS(-US)](https://ral.ucar.edu/solutions/products/camels).
+
 ### Download Daymet V4 dataset for basins in CAMELS
 
 Use hydrobench/app/download/download_daymet_camels_basin.py to download daymet grid data for the boundaries of basins in
diff --git a/docs/source/develop.rst b/docs/source/develop.rst
index 14d58e7..a7c4c90 100644
--- a/docs/source/develop.rst
+++ b/docs/source/develop.rst
@@ -20,7 +20,10 @@ Now the directories in "hydrobench" are:
 
 * app: scripts for downloading and preprocessing source data
 * data: this directory contains the data-source classes for source data;
-* daymet4basins: this is the directory for downloading daymet data for basins' boundaries;
+* daymet4basins: this is the directory for downloading and processing daymet data for camels-basins' boundaries;
+* ecmwf4basins: this is the directory for downloading and processing ERA5 data;
+* modis4basins: this is the directory for downloading and processing MODIS products data;
+* nldas4basins: this is the directory for downloading and processing NLDAS data;
 * pet: all methods for calculating potential evapotranspiration are here;
 * utils: all utility functions.
 
@@ -28,3 +31,14 @@ Conventions
 ----------------------
 * Please write unittest code for all functions.
 * Please write comments by English, and note we use the the same format in `this example <https://sphinxcontrib-napoleon.readthedocs.io/en/latest/example_numpy.html#example-numpy>`_.
+
+Update documents
+----------------------
+After adding new features or modifying existed code, please update the document:
+
+.. code-block:: shell
+
+    cd docs
+    sphinx-apidoc -o source/ ../hydrobench
+    make html
+
diff --git a/docs/source/hydrobench.app.daymet4basins.rst b/docs/source/hydrobench.app.daymet4basins.rst
index 2690cf2..0f954ea 100644
--- a/docs/source/hydrobench.app.daymet4basins.rst
+++ b/docs/source/hydrobench.app.daymet4basins.rst
@@ -4,6 +4,14 @@ hydrobench.app.daymet4basins package
 Submodules
 ----------
 
+hydrobench.app.daymet4basins.calculate\_basin\_mean\_forcing\_include\_pet module
+---------------------------------------------------------------------------------
+
+.. automodule:: hydrobench.app.daymet4basins.calculate_basin_mean_forcing_include_pet
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 hydrobench.app.daymet4basins.calculate\_basin\_mean\_values module
 ------------------------------------------------------------------
 
@@ -20,6 +28,14 @@ hydrobench.app.daymet4basins.cut\_nc\_to\_bounds\_add\_pet module
    :undoc-members:
    :show-inheritance:
 
+hydrobench.app.daymet4basins.pbm\_p\_pe\_q\_basin\_area module
+--------------------------------------------------------------
+
+.. automodule:: hydrobench.app.daymet4basins.pbm_p_pe_q_basin_area
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 hydrobench.app.daymet4basins.regrid\_daymet\_nc module
 ------------------------------------------------------
 
diff --git a/docs/source/hydrobench.app.download.rst b/docs/source/hydrobench.app.download.rst
index 5f30522..52aec2b 100644
--- a/docs/source/hydrobench.app.download.rst
+++ b/docs/source/hydrobench.app.download.rst
@@ -28,6 +28,14 @@ hydrobench.app.download.download\_daymet\_nldi\_camels\_basin module
    :undoc-members:
    :show-inheritance:
 
+hydrobench.app.download.download\_nldas\_hourly module
+------------------------------------------------------
+
+.. automodule:: hydrobench.app.download.download_nldas_hourly
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 hydrobench.app.download.download\_nldi module
 ---------------------------------------------
 
diff --git a/docs/source/hydrobench.app.modis4basins.rst b/docs/source/hydrobench.app.modis4basins.rst
new file mode 100644
index 0000000..5871f58
--- /dev/null
+++ b/docs/source/hydrobench.app.modis4basins.rst
@@ -0,0 +1,21 @@
+hydrobench.app.modis4basins package
+===================================
+
+Submodules
+----------
+
+hydrobench.app.modis4basins.trans\_modis\_et\_to\_camels\_format module
+-----------------------------------------------------------------------
+
+.. automodule:: hydrobench.app.modis4basins.trans_modis_et_to_camels_format
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Module contents
+---------------
+
+.. automodule:: hydrobench.app.modis4basins
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/source/hydrobench.app.nldas4basins.rst b/docs/source/hydrobench.app.nldas4basins.rst
new file mode 100644
index 0000000..e6f0c37
--- /dev/null
+++ b/docs/source/hydrobench.app.nldas4basins.rst
@@ -0,0 +1,21 @@
+hydrobench.app.nldas4basins package
+===================================
+
+Submodules
+----------
+
+hydrobench.app.nldas4basins.trans\_nldas\_to\_camels\_format module
+-------------------------------------------------------------------
+
+.. automodule:: hydrobench.app.nldas4basins.trans_nldas_to_camels_format
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Module contents
+---------------
+
+.. automodule:: hydrobench.app.nldas4basins
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/source/hydrobench.app.preprocess.rst b/docs/source/hydrobench.app.preprocess.rst
deleted file mode 100644
index 24ca261..0000000
--- a/docs/source/hydrobench.app.preprocess.rst
+++ /dev/null
@@ -1,37 +0,0 @@
-hydrobench.app.preprocess package
-=================================
-
-Submodules
-----------
-
-hydrobench.app.preprocess.calculate\_basin\_mean\_values module
----------------------------------------------------------------
-
-.. automodule:: hydrobench.app.preprocess.calculate_basin_mean_values
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-hydrobench.app.preprocess.calculate\_pet module
------------------------------------------------
-
-.. automodule:: hydrobench.app.preprocess.calculate_pet
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-hydrobench.app.preprocess.cut\_nc\_to\_bounds module
-----------------------------------------------------
-
-.. automodule:: hydrobench.app.preprocess.cut_nc_to_bounds
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: hydrobench.app.preprocess
-   :members:
-   :undoc-members:
-   :show-inheritance:
diff --git a/docs/source/hydrobench.app.rst b/docs/source/hydrobench.app.rst
index 5329843..3e53e5f 100644
--- a/docs/source/hydrobench.app.rst
+++ b/docs/source/hydrobench.app.rst
@@ -9,6 +9,8 @@ Subpackages
 
    hydrobench.app.daymet4basins
    hydrobench.app.download
+   hydrobench.app.modis4basins
+   hydrobench.app.nldas4basins
 
 Submodules
 ----------
diff --git a/docs/source/hydrobench.data.rst b/docs/source/hydrobench.data.rst
index 8912c42..68db710 100644
--- a/docs/source/hydrobench.data.rst
+++ b/docs/source/hydrobench.data.rst
@@ -20,6 +20,14 @@ hydrobench.data.data\_camels module
    :undoc-members:
    :show-inheritance:
 
+hydrobench.data.data\_daymet4camels module
+------------------------------------------
+
+.. automodule:: hydrobench.data.data_daymet4camels
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 hydrobench.data.stat module
 ---------------------------
 
diff --git a/docs/source/hydrobench.ecmwf4basins.rst b/docs/source/hydrobench.ecmwf4basins.rst
new file mode 100644
index 0000000..979fdf2
--- /dev/null
+++ b/docs/source/hydrobench.ecmwf4basins.rst
@@ -0,0 +1,21 @@
+hydrobench.ecmwf4basins package
+===============================
+
+Submodules
+----------
+
+hydrobench.ecmwf4basins.download\_era5\_land module
+---------------------------------------------------
+
+.. automodule:: hydrobench.ecmwf4basins.download_era5_land
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Module contents
+---------------
+
+.. automodule:: hydrobench.ecmwf4basins
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/source/hydrobench.modis4basins.rst b/docs/source/hydrobench.modis4basins.rst
new file mode 100644
index 0000000..d329eb7
--- /dev/null
+++ b/docs/source/hydrobench.modis4basins.rst
@@ -0,0 +1,53 @@
+hydrobench.modis4basins package
+===============================
+
+Submodules
+----------
+
+hydrobench.modis4basins.DAACDataDownload module
+-----------------------------------------------
+
+.. automodule:: hydrobench.modis4basins.DAACDataDownload
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+hydrobench.modis4basins.EarthdataLoginSetup module
+--------------------------------------------------
+
+.. automodule:: hydrobench.modis4basins.EarthdataLoginSetup
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+hydrobench.modis4basins.basin\_mod16a2v105\_process module
+----------------------------------------------------------
+
+.. automodule:: hydrobench.modis4basins.basin_mod16a2v105_process
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+hydrobench.modis4basins.basin\_pmlv2\_process module
+----------------------------------------------------
+
+.. automodule:: hydrobench.modis4basins.basin_pmlv2_process
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+hydrobench.modis4basins.get\_modis module
+-----------------------------------------
+
+.. automodule:: hydrobench.modis4basins.get_modis
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Module contents
+---------------
+
+.. automodule:: hydrobench.modis4basins
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/source/hydrobench.nldas4basins.rst b/docs/source/hydrobench.nldas4basins.rst
index 45362fe..b6d241b 100644
--- a/docs/source/hydrobench.nldas4basins.rst
+++ b/docs/source/hydrobench.nldas4basins.rst
@@ -12,6 +12,14 @@ hydrobench.nldas4basins.basin\_nldas\_process module
    :undoc-members:
    :show-inheritance:
 
+hydrobench.nldas4basins.download\_nldas module
+----------------------------------------------
+
+.. automodule:: hydrobench.nldas4basins.download_nldas
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 Module contents
 ---------------
 
diff --git a/docs/source/hydrobench.rst b/docs/source/hydrobench.rst
index f5bffb3..5844016 100644
--- a/docs/source/hydrobench.rst
+++ b/docs/source/hydrobench.rst
@@ -10,6 +10,8 @@ Subpackages
    hydrobench.app
    hydrobench.data
    hydrobench.daymet4basins
+   hydrobench.ecmwf4basins
+   hydrobench.modis4basins
    hydrobench.nldas4basins
    hydrobench.pet
    hydrobench.utils
diff --git a/environment.yml b/environment.yml
index d5458e2..94b4bc1 100644
--- a/environment.yml
+++ b/environment.yml
@@ -9,6 +9,7 @@ dependencies:
   - pynhd=0.11.0
   - pygeohydro=0.11.0
   - pydaymet=0.11.0
+  - openpyxl
   - cdsapi
   - cfgrib
   - tqdm
diff --git a/hydrobench/data/data_camels.py b/hydrobench/data/data_camels.py
index ae193fd..b117d4d 100644
--- a/hydrobench/data/data_camels.py
+++ b/hydrobench/data/data_camels.py
@@ -1,5 +1,8 @@
 import collections
+import fnmatch
 import os
+from typing import Union
+import tarfile
 import pandas as pd
 import numpy as np
 from pandas.core.dtypes.common import is_string_dtype, is_numeric_dtype
@@ -7,94 +10,546 @@
 from hydrobench.data.data_base import DataSourceBase
 from hydrobench.data.stat import cal_fdc
 from hydrobench.utils import hydro_utils
-from hydrobench.utils.hydro_utils import download_one_zip
+from hydrobench.utils.hydro_utils import download_one_zip, unzip_nested_zip
+
+CAMELS_NO_DATASET_ERROR_LOG = "We cannot read this dataset now. Please check if you choose the correct dataset:\n" \
+                              " [\"AUS\", \"BR\", \"CA\", \"CL\", \"GB\", \"US\", \"YR\"]"
+
+
+def time_intersect_dynamic_data(obs: np.array, date: np.array, t_range: list):
+    """
+    chose data from obs in the t_range
+
+    Parameters
+    ----------
+    obs
+        a np array
+    date
+        all periods for obs
+    t_range
+        the time range we need, such as ["1990-01-01","2000-01-01"]
+
+    Returns
+    -------
+    np.array
+        the chosen data
+    """
+    t_lst = hydro_utils.t_range_days(t_range)
+    nt = t_lst.shape[0]
+    if len(obs) != nt:
+        out = np.full([nt], np.nan)
+        [c, ind1, ind2] = np.intersect1d(date, t_lst, return_indices=True)
+        out[ind2] = obs[ind1]
+    else:
+        out = obs
+    return out
 
 
 class Camels(DataSourceBase):
-    def __init__(self, data_path, download=False):
+    def __init__(self, data_path, download=False, region: str = "US"):
+        """
+        Initialization for CAMELS series dataset
+
+        Parameters
+        ----------
+        data_path
+            where we put the dataset
+        download
+            if true, download
+        region
+            the default is CAMELS(-US), since it's the first CAMELS dataset.
+            Others now include: AUS, BR, CL, GB, YR
+        """
         super().__init__(data_path)
+        region_lst = ["AUS", "BR", "CA", "CE", "CL", "GB", "US", "YR"]
+        assert region in region_lst
+        self.region = region
         self.data_source_description = self.set_data_source_describe()
         if download:
             self.download_data_source()
         self.camels_sites = self.read_site_info()
 
     def get_name(self):
-        return "CAMELS"
+        return "CAMELS_" + self.region
 
-    def set_data_source_describe(self):
+    def set_data_source_describe(self) -> collections.OrderedDict:
+        """
+        Introduce the files in the dataset and list their location in the file system
+
+        Returns
+        -------
+        collections.OrderedDict
+            the description for a CAMELS dataset
+        """
         camels_db = self.data_source_dir
-        # shp file of basins
-        camels_shp_file = os.path.join(camels_db, "basin_set_full_res", "HCDN_nhru_final_671.shp")
-        # config of flow data
-        flow_dir = os.path.join(camels_db, "basin_timeseries_v1p2_metForcing_obsFlow", "basin_dataset_public_v1p2",
-                                "usgs_streamflow")
-        # forcing
-        forcing_dir = os.path.join(camels_db, "basin_timeseries_v1p2_metForcing_obsFlow", "basin_dataset_public_v1p2",
-                                   "basin_mean_forcing")
-        forcing_types = ["daymet", "maurer", "nldas"]
-        # attr
-        attr_dir = os.path.join(camels_db, "camels_attributes_v2.0", "camels_attributes_v2.0")
-        gauge_id_file = os.path.join(attr_dir, 'camels_name.txt')
-
-        download_url_lst = [
-            "https://ral.ucar.edu/sites/default/files/public/product-tool/camels-catchment-attributes-and-meteorology-for-large-sample-studies-dataset-downloads/camels_attributes_v2.0.zip",
-            "https://ral.ucar.edu/sites/default/files/public/product-tool/camels-catchment-attributes-and-meteorology-for-large-sample-studies-dataset-downloads/basin_set_full_res.zip",
-            "https://ral.ucar.edu/sites/default/files/public/product-tool/camels-catchment-attributes-and-meteorology-for-large-sample-studies-dataset-downloads/basin_timeseries_v1p2_metForcing_obsFlow.zip"]
-
-        return collections.OrderedDict(CAMELS_DIR=camels_db, CAMELS_FLOW_DIR=flow_dir,
-                                       CAMELS_FORCING_DIR=forcing_dir, CAMELS_FORCING_TYPE=forcing_types,
-                                       CAMELS_ATTR_DIR=attr_dir, CAMELS_GAUGE_FILE=gauge_id_file,
-                                       CAMELS_BASINS_SHP_FILE=camels_shp_file, CAMELS_DOWNLOAD_URL_LST=download_url_lst)
-
-    def download_data_source(self):
+        if self.region == "US":
+            # shp file of basins
+            camels_shp_file = os.path.join(camels_db, "basin_set_full_res", "HCDN_nhru_final_671.shp")
+            # config of flow data
+            flow_dir = os.path.join(camels_db, "basin_timeseries_v1p2_metForcing_obsFlow", "basin_dataset_public_v1p2",
+                                    "usgs_streamflow")
+            # forcing
+            forcing_dir = os.path.join(camels_db, "basin_timeseries_v1p2_metForcing_obsFlow",
+                                       "basin_dataset_public_v1p2",
+                                       "basin_mean_forcing")
+            forcing_types = ["daymet", "maurer", "nldas"]
+            # attr
+            attr_dir = os.path.join(camels_db, "camels_attributes_v2.0", "camels_attributes_v2.0")
+            gauge_id_file = os.path.join(attr_dir, 'camels_name.txt')
+            attr_key_lst = ['topo', 'clim', 'hydro', 'vege', 'soil', 'geol']
+            download_url_lst = [
+                "https://ral.ucar.edu/sites/default/files/public/product-tool/camels-catchment-attributes-and-meteorology-for-large-sample-studies-dataset-downloads/camels_attributes_v2.0.zip",
+                "https://ral.ucar.edu/sites/default/files/public/product-tool/camels-catchment-attributes-and-meteorology-for-large-sample-studies-dataset-downloads/basin_set_full_res.zip",
+                "https://ral.ucar.edu/sites/default/files/public/product-tool/camels-catchment-attributes-and-meteorology-for-large-sample-studies-dataset-downloads/basin_timeseries_v1p2_metForcing_obsFlow.zip"]
+
+            return collections.OrderedDict(CAMELS_DIR=camels_db, CAMELS_FLOW_DIR=flow_dir,
+                                           CAMELS_FORCING_DIR=forcing_dir, CAMELS_FORCING_TYPE=forcing_types,
+                                           CAMELS_ATTR_DIR=attr_dir, CAMELS_ATTR_KEY_LST=attr_key_lst,
+                                           CAMELS_GAUGE_FILE=gauge_id_file,
+                                           CAMELS_BASINS_SHP_FILE=camels_shp_file,
+                                           CAMELS_DOWNLOAD_URL_LST=download_url_lst)
+        elif self.region == "AUS":
+            # id and name
+            gauge_id_file = os.path.join(camels_db, "01_id_name_metadata", "01_id_name_metadata",
+                                         "id_name_metadata.csv")
+            # shp file of basins
+            camels_shp_file = os.path.join(camels_db, "02_location_boundary_area", "02_location_boundary_area", "shp",
+                                           "CAMELS_AUS_BasinOutlets_adopted.shp")
+            # config of flow data
+            flow_dir = os.path.join(camels_db, "03_streamflow", "03_streamflow")
+            # attr
+            attr_dir = os.path.join(camels_db, "04_attributes", "04_attributes")
+            # forcing
+            forcing_dir = os.path.join(camels_db, "05_hydrometeorology", "05_hydrometeorology")
+
+            return collections.OrderedDict(CAMELS_DIR=camels_db, CAMELS_FLOW_DIR=flow_dir,
+                                           CAMELS_FORCING_DIR=forcing_dir,
+                                           CAMELS_ATTR_DIR=attr_dir, CAMELS_GAUGE_FILE=gauge_id_file,
+                                           CAMELS_BASINS_SHP_FILE=camels_shp_file)
+        elif self.region == "BR":
+            # attr
+            attr_dir = os.path.join(camels_db, "01_CAMELS_BR_attributes", "01_CAMELS_BR_attributes")
+            # we don't need the location attr file
+            attr_key_lst = ['climate', 'geology', 'human_intervention', 'hydrology', 'land_cover',
+                            'quality_check', 'soil', 'topography']
+            # id and name, there are two types stations in CAMELS_BR, and we only chose the 897-stations version
+            gauge_id_file = os.path.join(attr_dir, "camels_br_topography.txt")
+            # shp file of basins
+            camels_shp_file = os.path.join(camels_db, "14_CAMELS_BR_catchment_boundaries",
+                                           "14_CAMELS_BR_catchment_boundaries", "camels_br_catchments.shp")
+            # config of flow data
+            flow_dir_m3s = os.path.join(camels_db, "02_CAMELS_BR_streamflow_m3s", "02_CAMELS_BR_streamflow_m3s")
+            flow_dir_mm_selected_catchments = os.path.join(camels_db, "03_CAMELS_BR_streamflow_mm_selected_catchments",
+                                                           "03_CAMELS_BR_streamflow_mm_selected_catchments")
+            flow_dir_simulated = os.path.join(camels_db, "04_CAMELS_BR_streamflow_simulated",
+                                              "04_CAMELS_BR_streamflow_simulated")
+
+            # forcing
+            forcing_dir_precipitation_chirps = os.path.join(camels_db, "05_CAMELS_BR_precipitation_chirps",
+                                                            "05_CAMELS_BR_precipitation_chirps")
+            forcing_dir_precipitation_mswep = os.path.join(camels_db, "06_CAMELS_BR_precipitation_mswep",
+                                                           "06_CAMELS_BR_precipitation_mswep")
+            forcing_dir_precipitation_cpc = os.path.join(camels_db, "07_CAMELS_BR_precipitation_cpc",
+                                                         "07_CAMELS_BR_precipitation_cpc")
+            forcing_dir_evapotransp_gleam = os.path.join(camels_db, "08_CAMELS_BR_evapotransp_gleam",
+                                                         "08_CAMELS_BR_evapotransp_gleam")
+            forcing_dir_evapotransp_mgb = os.path.join(camels_db, "09_CAMELS_BR_evapotransp_mgb",
+                                                       "09_CAMELS_BR_evapotransp_mgb")
+            forcing_dir_potential_evapotransp_gleam = os.path.join(camels_db,
+                                                                   "10_CAMELS_BR_potential_evapotransp_gleam",
+                                                                   "10_CAMELS_BR_potential_evapotransp_gleam")
+            forcing_dir_temperature_min_cpc = os.path.join(camels_db, "11_CAMELS_BR_temperature_min_cpc",
+                                                           "11_CAMELS_BR_temperature_min_cpc")
+            forcing_dir_temperature_mean_cpc = os.path.join(camels_db, "12_CAMELS_BR_temperature_mean_cpc",
+                                                            "12_CAMELS_BR_temperature_mean_cpc")
+            forcing_dir_temperature_max_cpc = os.path.join(camels_db, "13_CAMELS_BR_temperature_max_cpc",
+                                                           "13_CAMELS_BR_temperature_max_cpc")
+            return collections.OrderedDict(CAMELS_DIR=camels_db,
+                                           CAMELS_FLOW_DIR=[flow_dir_m3s, flow_dir_mm_selected_catchments,
+                                                            flow_dir_simulated],
+                                           CAMELS_FORCING_DIR=[forcing_dir_precipitation_chirps,
+                                                               forcing_dir_precipitation_mswep,
+                                                               forcing_dir_precipitation_cpc,
+                                                               forcing_dir_evapotransp_gleam,
+                                                               forcing_dir_evapotransp_mgb,
+                                                               forcing_dir_potential_evapotransp_gleam,
+                                                               forcing_dir_temperature_min_cpc,
+                                                               forcing_dir_temperature_mean_cpc,
+                                                               forcing_dir_temperature_max_cpc],
+                                           CAMELS_ATTR_DIR=attr_dir, CAMELS_ATTR_KEY_LST=attr_key_lst,
+                                           CAMELS_GAUGE_FILE=gauge_id_file,
+                                           CAMELS_BASINS_SHP_FILE=camels_shp_file)
+        elif self.region == "CL":
+            # attr
+            attr_dir = os.path.join(camels_db, "1_CAMELScl_attributes")
+            attr_file = os.path.join(attr_dir, "1_CAMELScl_attributes.txt")
+            # shp file of basins
+            camels_shp_file = os.path.join(camels_db, "CAMELScl_catchment_boundaries", "catchments_camels_cl_v1.3.shp")
+            # config of flow data
+            flow_dir_m3s = os.path.join(camels_db, "2_CAMELScl_streamflow_m3s")
+            flow_dir_mm = os.path.join(camels_db, "3_CAMELScl_streamflow_mm")
+
+            # forcing
+            forcing_dir_precip_cr2met = os.path.join(camels_db, "4_CAMELScl_precip_cr2met")
+            forcing_dir_precip_chirps = os.path.join(camels_db, "5_CAMELScl_precip_chirps")
+            forcing_dir_precip_mswep = os.path.join(camels_db, "6_CAMELScl_precip_mswep")
+            forcing_dir_precip_tmpa = os.path.join(camels_db, "7_CAMELScl_precip_tmpa")
+            forcing_dir_tmin_cr2met = os.path.join(camels_db, "8_CAMELScl_tmin_cr2met")
+            forcing_dir_tmax_cr2met = os.path.join(camels_db, "9_CAMELScl_tmax_cr2met")
+            forcing_dir_tmean_cr2met = os.path.join(camels_db, "10_CAMELScl_tmean_cr2met")
+            forcing_dir_pet_8d_modis = os.path.join(camels_db, "11_CAMELScl_pet_8d_modis")
+            forcing_dir_pet_hargreaves = os.path.join(camels_db, "12_CAMELScl_pet_hargreaves", )
+            forcing_dir_swe = os.path.join(camels_db, "13_CAMELScl_swe")
+            return collections.OrderedDict(CAMELS_DIR=camels_db,
+                                           CAMELS_FLOW_DIR=[flow_dir_m3s, flow_dir_mm],
+                                           CAMELS_FORCING_DIR=[forcing_dir_precip_cr2met, forcing_dir_precip_chirps,
+                                                               forcing_dir_precip_mswep, forcing_dir_precip_tmpa,
+                                                               forcing_dir_tmin_cr2met, forcing_dir_tmax_cr2met,
+                                                               forcing_dir_tmean_cr2met, forcing_dir_pet_8d_modis,
+                                                               forcing_dir_pet_hargreaves, forcing_dir_swe],
+                                           CAMELS_ATTR_DIR=attr_dir,
+                                           CAMELS_GAUGE_FILE=attr_file,
+                                           CAMELS_BASINS_SHP_FILE=camels_shp_file)
+        elif self.region == "GB":
+            # shp file of basins
+            camels_shp_file = os.path.join(camels_db, "8344e4f3-d2ea-44f5-8afa-86d2987543a9",
+                                           "8344e4f3-d2ea-44f5-8afa-86d2987543a9", "data",
+                                           "CAMELS_GB_catchment_boundaries", "CAMELS_GB_catchment_boundaries.shp")
+            # flow and forcing data are in a same file
+            flow_dir = os.path.join(camels_db, "8344e4f3-d2ea-44f5-8afa-86d2987543a9",
+                                    "8344e4f3-d2ea-44f5-8afa-86d2987543a9", "data", "timeseries")
+            forcing_dir = flow_dir
+            # attr
+            attr_dir = os.path.join(camels_db, "8344e4f3-d2ea-44f5-8afa-86d2987543a9",
+                                    "8344e4f3-d2ea-44f5-8afa-86d2987543a9", "data")
+            gauge_id_file = os.path.join(attr_dir, 'CAMELS_GB_hydrometry_attributes.csv')
+            attr_key_lst = ["climatic", "humaninfluence", "hydrogeology", "hydrologic", "hydrometry", "landcover",
+                            "soil", "topographic"]
+
+            return collections.OrderedDict(CAMELS_DIR=camels_db, CAMELS_FLOW_DIR=flow_dir,
+                                           CAMELS_FORCING_DIR=forcing_dir, CAMELS_ATTR_DIR=attr_dir,
+                                           CAMELS_ATTR_KEY_LST=attr_key_lst, CAMELS_GAUGE_FILE=gauge_id_file,
+                                           CAMELS_BASINS_SHP_FILE=camels_shp_file)
+        elif self.region == "YR":
+            # shp files of basins
+            camels_shp_files_dir = os.path.join(camels_db, "9_Normal_Camels_YR", "Normal_Camels_YR_basin_boundary")
+            # attr, flow and forcing data are all in the same dir. each basin has one dir.
+            flow_dir = os.path.join(camels_db, "9_Normal_Camels_YR", "1_Normal_Camels_YR_basin_data")
+            forcing_dir = flow_dir
+            attr_dir = flow_dir
+            # no gauge id file for CAMELS_YR; natural_watersheds.txt showed unregulated basins in CAMELS_YR
+            gauge_id_file = os.path.join(camels_db, "9_Normal_Camels_YR", "natural_watersheds.txt")
+            return collections.OrderedDict(CAMELS_DIR=camels_db, CAMELS_FLOW_DIR=flow_dir,
+                                           CAMELS_FORCING_DIR=forcing_dir, CAMELS_ATTR_DIR=attr_dir,
+                                           CAMELS_GAUGE_FILE=gauge_id_file, CAMELS_BASINS_SHP_DIR=camels_shp_files_dir)
+        elif self.region == "CA":
+            # shp file of basins
+            camels_shp_files_dir = os.path.join(camels_db, "CANOPEX_BOUNDARIES")
+            # config of flow data
+            flow_dir = os.path.join(camels_db, "CANOPEX_NRCAN_ASCII", "CANOPEX_NRCAN_ASCII")
+            forcing_dir = flow_dir
+            # There is no attr data in CANOPEX, hence we use attr from HYSET -- https://osf.io/7fn4c/
+            attr_dir = camels_db
+
+            gauge_id_file = os.path.join(camels_db, 'STATION_METADATA.xlsx')
+
+            return collections.OrderedDict(CAMELS_DIR=camels_db, CAMELS_FLOW_DIR=flow_dir,
+                                           CAMELS_FORCING_DIR=forcing_dir,
+                                           CAMELS_ATTR_DIR=attr_dir,
+                                           CAMELS_GAUGE_FILE=gauge_id_file,
+                                           CAMELS_BASINS_SHP_DIR=camels_shp_files_dir)
+        elif self.region == "CE":
+            # We use A_basins_total_upstrm
+            # shp file of basins
+            camels_shp_file = os.path.join(camels_db, "2_LamaH-CE_daily", "A_basins_total_upstrm", "3_shapefiles",
+                                           "Basins_A.shp")
+            # config of flow data
+            flow_dir = os.path.join(camels_db, "2_LamaH-CE_daily", "D_gauges", "2_timeseries", "daily")
+            forcing_dir = os.path.join(camels_db, "2_LamaH-CE_daily", "A_basins_total_upstrm", "2_timeseries", "daily")
+            attr_dir = os.path.join(camels_db, "2_LamaH-CE_daily", "A_basins_total_upstrm", "1_attributes")
+
+            gauge_id_file = os.path.join(camels_db, "2_LamaH-CE_daily", "D_gauges", "1_attributes",
+                                         "Gauge_attributes.csv")
+
+            return collections.OrderedDict(CAMELS_DIR=camels_db, CAMELS_FLOW_DIR=flow_dir,
+                                           CAMELS_FORCING_DIR=forcing_dir,
+                                           CAMELS_ATTR_DIR=attr_dir,
+                                           CAMELS_GAUGE_FILE=gauge_id_file,
+                                           CAMELS_BASINS_SHP_FILE=camels_shp_file)
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
+
+    def download_data_source(self) -> None:
+        """
+        Download CAMELS dataset.
+
+        Now we only support CAMELS-US's downloading.
+        For others, please download it manually and put all files of a CAMELS dataset in one directory.
+        For example, all files of CAMELS_AUS should be put in "camels_aus" directory
+
+        Returns
+        -------
+        None
+        """
         camels_config = self.data_source_description
-        if not os.path.isdir(camels_config["CAMELS_DIR"]):
-            os.makedirs(camels_config["CAMELS_DIR"])
-        [download_one_zip(attr_url, camels_config["CAMELS_DIR"]) for attr_url in
-         camels_config["CAMELS_DOWNLOAD_URL_LST"] if
-         not os.path.isfile(os.path.join(camels_config["CAMELS_DIR"], attr_url.split("/")[-1]))]
-        print("The CAMELS data have been downloaded!")
+        if self.region == "US":
+            if not os.path.isdir(camels_config["CAMELS_DIR"]):
+                os.makedirs(camels_config["CAMELS_DIR"])
+            [download_one_zip(attr_url, camels_config["CAMELS_DIR"]) for attr_url in
+             camels_config["CAMELS_DOWNLOAD_URL_LST"] if
+             not os.path.isfile(os.path.join(camels_config["CAMELS_DIR"], attr_url.split("/")[-1]))]
+            print("The CAMELS_US data have been downloaded!")
+        print("Please download it manually and put all files of a CAMELS dataset in the CAMELS_DIR directory.")
+        print("We unzip all files now.")
+        if self.region == "CE":
+            # We only use CE's dauly files now and it is tar.gz formatting
+            file = tarfile.open(os.path.join(camels_config["CAMELS_DIR"], "2_LamaH-CE_daily.tar.gz"))
+            # extracting file
+            file.extractall(os.path.join(camels_config["CAMELS_DIR"], "2_LamaH-CE_daily"))
+            file.close()
+        for f_name in os.listdir(camels_config["CAMELS_DIR"]):
+            if fnmatch.fnmatch(f_name, '*.zip'):
+                unzip_dir = os.path.join(camels_config["CAMELS_DIR"], f_name[0:-4])
+                file_name = os.path.join(camels_config["CAMELS_DIR"], f_name)
+                unzip_nested_zip(file_name, unzip_dir)
+
+    def read_site_info(self) -> pd.DataFrame:
+        """
+        Read the basic information of gages in a CAMELS dataset
+
+        Returns
+        -------
+        pd.DataFrame
+            basic info of gages
+        """
+        camels_file = self.data_source_description["CAMELS_GAUGE_FILE"]
+        if self.region == "US":
+            data = pd.read_csv(camels_file, sep=';', dtype={"gauge_id": str, "huc_02": str})
+        elif self.region == "AUS":
+            data = pd.read_csv(camels_file, sep=',', dtype={"station_id": str})
+        elif self.region == "BR":
+            data = pd.read_csv(camels_file, sep='\s+', dtype={"gauge_id": str})
+        elif self.region == "CL":
+            data = pd.read_csv(camels_file, sep='\t', index_col=0)
+        elif self.region == "GB":
+            data = pd.read_csv(camels_file, sep=',', dtype={"gauge_id": str})
+        elif self.region == "YR":
+            dirs_ = os.listdir(self.data_source_description["CAMELS_ATTR_DIR"])
+            data = pd.DataFrame({"gauge_id": dirs_})
+        elif self.region == "CA":
+            data = pd.read_excel(camels_file)
+        elif self.region == "CE":
+            data = pd.read_csv(camels_file, sep=';')
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
+        return data
 
     def get_constant_cols(self) -> np.array:
-        """all readable attrs in CAMELS"""
+        """
+        all readable attrs in CAMELS
+
+        Returns
+        -------
+        np.array
+            attribute types
+        """
         data_folder = self.data_source_description["CAMELS_ATTR_DIR"]
-        var_dict = dict()
-        var_lst = list()
-        key_lst = ['topo', 'clim', 'hydro', 'vege', 'soil', 'geol']
-        for key in key_lst:
-            data_file = os.path.join(data_folder, 'camels_' + key + '.txt')
-            data_temp = pd.read_csv(data_file, sep=';')
-            var_lst_temp = list(data_temp.columns[1:])
-            var_dict[key] = var_lst_temp
-            var_lst.extend(var_lst_temp)
-        return np.array(var_lst)
+        if self.region == "US":
+            var_dict = dict()
+            var_lst = list()
+            key_lst = self.data_source_description["CAMELS_ATTR_KEY_LST"]
+            for key in key_lst:
+                data_file = os.path.join(data_folder, 'camels_' + key + '.txt')
+                data_temp = pd.read_csv(data_file, sep=';')
+                var_lst_temp = list(data_temp.columns[1:])
+                var_dict[key] = var_lst_temp
+                var_lst.extend(var_lst_temp)
+            return np.array(var_lst)
+        elif self.region == "AUS":
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_DIR"],
+                                         "CAMELS_AUS_Attributes-Indices_MasterTable.csv")
+            camels_aus_attr_indices_data = pd.read_csv(attr_all_file, sep=',')
+            # exclude station id
+            return camels_aus_attr_indices_data.columns.values[1:]
+        elif self.region == "BR":
+            var_dict = dict()
+            var_lst = list()
+            key_lst = self.data_source_description["CAMELS_ATTR_KEY_LST"]
+            for key in key_lst:
+                data_file = os.path.join(data_folder, 'camels_br_' + key + '.txt')
+                data_temp = pd.read_csv(data_file, sep='\s+')
+                var_lst_temp = list(data_temp.columns[1:])
+                var_dict[key] = var_lst_temp
+                var_lst.extend(var_lst_temp)
+            return np.array(var_lst)
+        elif self.region == "CL":
+            camels_cl_attr_data = self.camels_sites
+            # exclude station id
+            return camels_cl_attr_data.index.values
+        elif self.region == "GB":
+            var_dict = dict()
+            var_lst = list()
+            key_lst = self.data_source_description["CAMELS_ATTR_KEY_LST"]
+            for key in key_lst:
+                data_file = os.path.join(data_folder, "CAMELS_GB_" + key + '_attributes.csv')
+                data_temp = pd.read_csv(data_file, sep=',')
+                var_lst_temp = list(data_temp.columns[1:])
+                var_dict[key] = var_lst_temp
+                var_lst.extend(var_lst_temp)
+            return np.array(var_lst)
+        elif self.region == "YR":
+            attr_json_file = os.path.join(self.data_source_description["CAMELS_ATTR_DIR"], "0000", "attributes.json")
+            attr_json = hydro_utils.unserialize_json_ordered(attr_json_file)
+            return np.array(list(attr_json.keys()))
+        elif self.region == "CA":
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_DIR"], "HYSETS_watershed_properties.txt")
+            canopex_attr_indices_data = pd.read_csv(attr_all_file, sep=';')
+            # exclude HYSETS watershed id
+            return canopex_attr_indices_data.columns.values[1:]
+        elif self.region == "CE":
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_ATTR_DIR"], "Catchment_attributes.csv")
+            lamah_ce_attr_indices_data = pd.read_csv(attr_all_file, sep=';')
+            return lamah_ce_attr_indices_data.columns.values[1:]
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
 
     def get_relevant_cols(self) -> np.array:
-        # TODO: now only these 7 forcings
-        return np.array(['dayl', 'prcp', 'srad', 'swe', 'tmax', 'tmin', 'vp'])
+        """
+        all readable forcing types
+
+        Returns
+        -------
+        np.array
+            forcing types
+        """
+        if self.region == "US":
+            return np.array(["dayl", "prcp", "srad", "swe", "tmax", "tmin", "vp"])
+        elif self.region == "AUS":
+            forcing_types = []
+            for root, dirs, files in os.walk(self.data_source_description["CAMELS_FORCING_DIR"]):
+                if root == self.data_source_description["CAMELS_FORCING_DIR"]:
+                    continue
+                for file in files:
+                    forcing_types.append(file[:-4])
+            return np.array(forcing_types)
+        elif self.region == "BR":
+            return np.array([forcing_dir.split(os.sep)[-1][13:] for forcing_dir in
+                             self.data_source_description["CAMELS_FORCING_DIR"]])
+        elif self.region == "CL":
+            return np.array(["_".join(forcing_dir.split(os.sep)[-1].split("_")[2:]) for forcing_dir in
+                             self.data_source_description["CAMELS_FORCING_DIR"]])
+        elif self.region == "GB":
+            return np.array(["precipitation", "pet", "temperature", "peti", "humidity", "shortwave_rad", "longwave_rad",
+                             "windspeed"])
+        elif self.region == "YR":
+            return np.array(
+                ["pre", "evp", "gst_mean", "prs_mean", "tem_mean", "rhu", "win_mean", "gst_min", "prs_min", "tem_min",
+                 "gst_max", "prs_max", "tem_max", "ssd", "win_max"])
+        elif self.region == "CA":
+            # Although there is climatic potential evaporation item, CANOPEX does not have any PET data
+            return np.array(["prcp", "tmax", "tmin"])
+        elif self.region == "CE":
+            # Although there is climatic potential evaporation item, CANOPEX does not have any PET data
+            return np.array(["2m_temp_max", "2m_temp_mean", "2m_temp_min", "2m_dp_temp_max", "2m_dp_temp_mean",
+                             "2m_dp_temp_min", "10m_wind_u", "10m_wind_v", "fcst_alb", "lai_high_veg", "lai_low_veg",
+                             "swe", "surf_net_solar_rad_max", "surf_net_solar_rad_mean", "surf_net_therm_rad_max",
+                             "surf_net_therm_rad_mean", "surf_press", "total_et", "prec", "volsw_123", "volsw_4"])
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
 
     def get_target_cols(self) -> np.array:
-        # TODO: now only usgsFlow
-        return np.array(["usgsFlow"])
+        """
+        For CAMELS, the target vars are streamflows
+
+        Returns
+        -------
+        np.array
+            streamflow types
+        """
+        if self.region == "US":
+            return np.array(["usgsFlow"])
+        elif self.region == "AUS":
+            # QualityCodes are not streamflow data.
+            # MLd means "1 Megaliters Per Day"; 1 MLd = 0.011574074074074 cubic-meters-per-second
+            # mmd means "mm/day"
+            return np.array(
+                ["streamflow_MLd", "streamflow_MLd_inclInfilled", "streamflow_mmd", "streamflow_QualityCodes"])
+        elif self.region == "BR":
+            return np.array(
+                [flow_dir.split(os.sep)[-1][13:] for flow_dir in self.data_source_description["CAMELS_FLOW_DIR"]])
+        elif self.region == "CL":
+            return np.array(
+                [flow_dir.split(os.sep)[-1][11:] for flow_dir in self.data_source_description["CAMELS_FLOW_DIR"]])
+        elif self.region == "GB":
+            return np.array(["discharge_spec", "discharge_vol"])
+        elif self.region == "YR":
+            return np.array(["normalized_q"])
+        elif self.region == "CA":
+            return np.array(["discharge"])
+        elif self.region == "CE":
+            return np.array(["qobs"])
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
 
     def get_other_cols(self) -> dict:
         return {"FDC": {"time_range": ["1980-01-01", "2000-01-01"], "quantile_num": 100}}
 
-    def read_site_info(self) -> pd.DataFrame:
-        camels_file = self.data_source_description["CAMELS_GAUGE_FILE"]
-        data = pd.read_csv(camels_file, sep=';', dtype={"gauge_id": str, "huc_02": str})
-        return data
-
-    def read_object_ids(self, object_params=None) -> np.array:
-        return self.camels_sites["gauge_id"].values
+    def read_object_ids(self, **kwargs) -> np.array:
+        """
+        read station ids
 
-    def read_basin_area(self, object_ids) -> np.array:
-        return self.read_constant_cols(object_ids, ['area_gages2'], is_return_dict=False)
+        Parameters
+        ----------
+        **kwargs
+            optional params if needed
 
-    def read_mean_prep(self, object_ids) -> np.array:
-        return self.read_constant_cols(object_ids, ['p_mean'], is_return_dict=False)
+        Returns
+        -------
+        np.array
+            gage/station ids
+        """
+        if self.region in ["BR", "GB", "US", "YR"]:
+            return self.camels_sites["gauge_id"].values
+        elif self.region == "AUS":
+            return self.camels_sites["station_id"].values
+        elif self.region == "CL":
+            station_ids = self.camels_sites.columns.values
+            # for 7-digit id, replace the space with 0 to get a 8-digit id
+            cl_station_ids = [station_id.split(" ")[-1].zfill(8) for station_id in station_ids]
+            return np.array(cl_station_ids)
+        elif self.region == "CA":
+            ids = self.camels_sites["STATION_ID"].values
+            id_strs = [id_.split("'")[1] for id_ in ids]
+            # although there are 698 sites, there are only 611 sites with attributes data.
+            # Hence we only use 611 sites now
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_DIR"], "HYSETS_watershed_properties.txt")
+            canopex_attr_data = pd.read_csv(attr_all_file, sep=';')
+            return np.intersect1d(id_strs, canopex_attr_data["Official_ID"].values)
+        elif self.region == "CE":
+            # Not all basins have attributes, so we just chose those with attrs
+            ids = self.camels_sites["ID"].values
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_ATTR_DIR"], "Catchment_attributes.csv")
+            attr_data = pd.read_csv(attr_all_file, sep=';')
+            return np.intersect1d(ids, attr_data["ID"].values)
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
 
     def read_usgs_gage(self, usgs_id, t_range):
+        """
+        read streamflow data of a station from CAMELS-US
+
+        Parameters
+        ----------
+        usgs_id
+            the station id
+        t_range
+            the time range, for example, ["1990-01-01", "2000-01-01"]
+
+        Returns
+        -------
+        np.array
+            streamflow data of one station for a given time range
+        """
         print("reading %s streamflow data", usgs_id)
         gage_id_df = self.camels_sites
         huc = gage_id_df[gage_id_df["gauge_id"] == usgs_id]["huc_02"].values[0]
@@ -115,13 +570,187 @@ def read_usgs_gage(self, usgs_id, t_range):
             out = obs
         return out
 
-    def read_target_cols(self, usgs_id_lst=None, t_range=None, target_cols=None, **kwargs):
-        # only one output for camels: streamflow, just read it
-        nt = hydro_utils.t_range_days(t_range).shape[0]
-        y = np.empty([len(usgs_id_lst), nt])
-        for k in range(len(usgs_id_lst)):
-            data_obs = self.read_usgs_gage(usgs_id_lst[k], t_range)
-            y[k, :] = data_obs
+    def read_br_gage_flow(self, gage_id, t_range, flow_type):
+        """
+        Read gage's streamflow from CAMELS-BR
+
+        Parameters
+        ----------
+        gage_id
+            the station id
+        t_range
+            the time range, for example, ["1990-01-01", "2000-01-01"]
+        flow_type
+            "streamflow_m3s" or "streamflow_mm_selected_catchments" or "streamflow_simulated"
+
+        Returns
+        -------
+        np.array
+            streamflow data of one station for a given time range
+        """
+        dir_ = [flow_dir for flow_dir in self.data_source_description["CAMELS_FLOW_DIR"] if flow_type in flow_dir][0]
+        if flow_type == "streamflow_mm_selected_catchments":
+            flow_type = "streamflow_mm"
+        elif flow_type == "streamflow_simulated":
+            flow_type = "simulated_streamflow"
+        gage_file = os.path.join(dir_, gage_id + "_" + flow_type + ".txt")
+        data_temp = pd.read_csv(gage_file, sep=r'\s+')
+        obs = data_temp.iloc[:, 3].values
+        obs[obs < 0] = np.nan
+        df_date = data_temp[['year', 'month', 'day']]
+        date = pd.to_datetime(df_date).values.astype('datetime64[D]')
+        out = time_intersect_dynamic_data(obs, date, t_range)
+        return out
+
+    def read_gb_gage_flow_forcing(self, gage_id, t_range, var_type):
+        """
+        Read gage's streamflow or forcing from CAMELS-GB
+
+        Parameters
+        ----------
+        gage_id
+            the station id
+        t_range
+            the time range, for example, ["1990-01-01", "2000-01-01"]
+        var_type
+            flow type: "discharge_spec" or "discharge_vol"
+            forcing type: "precipitation", "pet", "temperature", "peti", "humidity", "shortwave_rad", "longwave_rad",
+                          "windspeed"
+
+        Returns
+        -------
+        np.array
+            streamflow or forcing data of one station for a given time range
+        """
+        gage_file = os.path.join(self.data_source_description["CAMELS_FLOW_DIR"],
+                                 "CAMELS_GB_hydromet_timeseries_" + gage_id + "_19701001-20150930.csv")
+        data_temp = pd.read_csv(gage_file, sep=",")
+        obs = data_temp[var_type].values
+        if var_type in ["discharge_spec", "discharge_vol"]:
+            obs[obs < 0] = np.nan
+        date = pd.to_datetime(data_temp["date"]).values.astype('datetime64[D]')
+        out = time_intersect_dynamic_data(obs, date, t_range)
+        return out
+
+    def read_target_cols(self, gage_id_lst: Union[list, np.array] = None, t_range: list = None,
+                         target_cols: Union[list, np.array] = None, **kwargs) -> np.array:
+        """
+        read target values; for CAMELS, they are streamflows
+        default target_cols is an one-value list
+
+        Parameters
+        ----------
+        gage_id_lst
+            station ids
+        t_range
+            the time range, for example, ["1990-01-01", "2000-01-01"]
+        target_cols
+            the default is None, but we neea at least one default target.
+            For CAMELS-US, it is ["usgsFlow"];
+            for CAMELS-AUS, it's ["streamflow_mmd"]
+            for CAMELS-AUS, it's ["streamflow_m3s"]
+        kwargs
+            some other params if needed
+
+        Returns
+        -------
+        np.array
+            streamflow data, 3-dim [station, time, streamflow]
+        """
+        if target_cols is None:
+            return np.array([])
+        else:
+            nf = len(target_cols)
+        t_range_list = hydro_utils.t_range_days(t_range)
+        nt = t_range_list.shape[0]
+        y = np.empty([len(gage_id_lst), nt, nf])
+        if self.region == "US":
+            for k in range(len(gage_id_lst)):
+                data_obs = self.read_usgs_gage(gage_id_lst[k], t_range)
+                # For CAMELS-US, only ["usgsFlow"]
+                y[k, :, 0] = data_obs
+        elif self.region == "AUS":
+            for k in range(len(target_cols)):
+                flow_data = pd.read_csv(
+                    os.path.join(self.data_source_description["CAMELS_FLOW_DIR"], target_cols[k] + ".csv"))
+                df_date = flow_data[['year', 'month', 'day']]
+                date = pd.to_datetime(df_date).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                chosen_data = flow_data[gage_id_lst].values[ind1, :]
+                chosen_data[chosen_data < 0] = np.nan
+                y[:, ind2, k] = chosen_data.T
+        elif self.region == "BR":
+            for j in range(len(target_cols)):
+                for k in range(len(gage_id_lst)):
+                    data_obs = self.read_br_gage_flow(gage_id_lst[k], t_range, target_cols[j])
+                    y[k, :, j] = data_obs
+        elif self.region == "CL":
+            for k in range(len(target_cols)):
+                if target_cols[k] == "streamflow_m3s":
+                    flow_data = pd.read_csv(os.path.join(self.data_source_description["CAMELS_FLOW_DIR"][0],
+                                                         "2_CAMELScl_streamflow_m3s.txt"), sep="\t", index_col=0)
+                elif target_cols[k] == "streamflow_mm":
+                    flow_data = pd.read_csv(os.path.join(self.data_source_description["CAMELS_FLOW_DIR"][1],
+                                                         "3_CAMELScl_streamflow_mm.txt"), sep="\t", index_col=0)
+                else:
+                    raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
+                date = pd.to_datetime(flow_data.index.values).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                station_ids = self.read_object_ids()
+                assert (all(x < y for x, y in zip(station_ids, station_ids[1:])))
+                [s, ind3, ind4] = np.intersect1d(station_ids, gage_id_lst, return_indices=True)
+                chosen_data = flow_data.iloc[ind1, ind3].replace("\s+", np.nan, regex=True)
+                chosen_data = chosen_data.astype(float)
+                chosen_data[chosen_data < 0] = np.nan
+                y[:, ind2, k] = chosen_data.values.T
+        elif self.region == "GB":
+            for j in range(len(target_cols)):
+                for k in range(len(gage_id_lst)):
+                    data_obs = self.read_gb_gage_flow_forcing(gage_id_lst[k], t_range, target_cols[j])
+                    y[k, :, j] = data_obs
+        elif self.region == "YR":
+            for k in range(len(gage_id_lst)):
+                # only one streamflow type: normalized_q
+                flow_file = os.path.join(self.data_source_description["CAMELS_FLOW_DIR"], gage_id_lst[k],
+                                         target_cols[0] + ".csv")
+                flow_data = pd.read_csv(flow_file, sep=",")
+                date = pd.to_datetime(flow_data["date"]).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                # flow data has been normalized, so we don't set negative values NaN
+                y[k, ind2, 0] = flow_data["q"].values[ind1]
+        elif self.region == "CA":
+            for k in range(len(gage_id_lst)):
+                # only one streamflow type: discharge
+                canopex_id = self.camels_sites[
+                    self.camels_sites["STATION_ID"] == "'" + gage_id_lst[k] + "'"]["CANOPEX_ID"].values[0]
+                flow_file = os.path.join(self.data_source_description["CAMELS_FLOW_DIR"], str(canopex_id) + ".dly")
+                read_flow_file = pd.read_csv(flow_file, header=None).values.tolist()
+                flow_data = []
+                flow_date = []
+                for one_site in read_flow_file:
+                    flow_date.append(hydro_utils.t2dt(int(one_site[0][:8].replace(" ", "0"))))
+                    all_data = one_site[0].split(" ")
+                    real_data = [one_data for one_data in all_data if one_data != ""]
+                    flow_data.append(float(real_data[-3]))
+                date = pd.to_datetime(flow_date).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                obs = np.array(flow_data)
+                obs[obs < 0] = np.nan
+                y[k, ind2, 0] = obs[ind1]
+        elif self.region == "CE":
+            for k in range(len(gage_id_lst)):
+                flow_file = os.path.join(self.data_source_description["CAMELS_FLOW_DIR"],
+                                         "ID_" + str(gage_id_lst[k]) + ".csv")
+                flow_data = pd.read_csv(flow_file, sep=";")
+                df_date = flow_data[["YYYY", "MM", "DD"]]
+                df_date.columns = ['year', 'month', 'day']
+                date = pd.to_datetime(df_date).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                obs = flow_data["qobs"].values
+                obs[obs < 0] = np.nan
+                y[k, ind2, 0] = obs[ind1]
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
         return y
 
     def read_forcing_gage(self, usgs_id, var_lst, t_range_list, forcing_type='daymet'):
@@ -149,29 +778,186 @@ def read_forcing_gage(self, usgs_id, var_lst, t_range_list, forcing_type='daymet
 
         for k in range(nf):
             ind = forcing_lst.index(var_lst[k])
-            out[:, k] = data_temp[ind].values[ind1]
+            out[ind2, k] = data_temp[ind].values[ind1]
         return out
 
-    def read_relevant_cols(self, usgs_id_lst=None, t_range=None, var_lst=None, forcing_type="daymet"):
+    def read_br_basin_forcing(self, gage_id, t_range, var_type) -> np.array:
+        """
+        Read one forcing data for a basin in CAMELS_BR
+
+        Parameters
+        ----------
+        gage_id
+            basin id
+        t_range
+            the time range, for example, ["1995-01-01", "2005-01-01"]
+        var_type
+            the forcing variable type
+
+        Returns
+        -------
+        np.array
+            one type forcing data of a basin in a given time range
+        """
+        dir_ = [_dir for _dir in self.data_source_description["CAMELS_FORCING_DIR"] if var_type in _dir][0]
+        if var_type in ['temperature_min_cpc', 'temperature_mean_cpc', 'temperature_max_cpc']:
+            var_type = var_type[:-4]
+        gage_file = os.path.join(dir_, gage_id + "_" + var_type + ".txt")
+        data_temp = pd.read_csv(gage_file, sep=r'\s+')
+        obs = data_temp.iloc[:, 3].values
+        df_date = data_temp[['year', 'month', 'day']]
+        date = pd.to_datetime(df_date).values.astype('datetime64[D]')
+        out = time_intersect_dynamic_data(obs, date, t_range)
+        return out
+
+    def read_relevant_cols(self, gage_id_lst: list = None, t_range: list = None, var_lst: list = None,
+                           forcing_type="daymet") -> np.array:
+        """
+        Read forcing data
+
+        Parameters
+        ----------
+        gage_id_lst
+            station ids
+        t_range
+            the time range, for example, ["1990-01-01", "2000-01-01"]
+        var_lst
+            forcing variable types
+        forcing_type
+            now only for CAMELS-US, there are three types: daymet, nldas, maurer
+        Returns
+        -------
+        np.array
+            forcing data
+        """
         t_range_list = hydro_utils.t_range_days(t_range)
         nt = t_range_list.shape[0]
-        x = np.empty([len(usgs_id_lst), nt, len(var_lst)])
-        for k in range(len(usgs_id_lst)):
-            data = self.read_forcing_gage(usgs_id_lst[k], var_lst, t_range_list, forcing_type=forcing_type)
-            x[k, :, :] = data
+        x = np.empty([len(gage_id_lst), nt, len(var_lst)])
+        if self.region == "US":
+            for k in range(len(gage_id_lst)):
+                data = self.read_forcing_gage(gage_id_lst[k], var_lst, t_range_list, forcing_type=forcing_type)
+                x[k, :, :] = data
+        elif self.region == "AUS":
+            for k in range(len(var_lst)):
+                if "precipitation_" in var_lst[k]:
+                    forcing_dir = os.path.join(self.data_source_description["CAMELS_FORCING_DIR"],
+                                               "01_precipitation_timeseries")
+                elif "et_" in var_lst[k] or "evap_" in var_lst[k]:
+                    forcing_dir = os.path.join(self.data_source_description["CAMELS_FORCING_DIR"],
+                                               "02_EvaporativeDemand_timeseries")
+                else:
+                    if "_AWAP" in var_lst[k]:
+                        forcing_dir = os.path.join(self.data_source_description["CAMELS_FORCING_DIR"], "03_Other",
+                                                   "AWAP")
+                    elif "_SILO" in var_lst[k]:
+                        forcing_dir = os.path.join(self.data_source_description["CAMELS_FORCING_DIR"], "03_Other",
+                                                   "SILO")
+                    else:
+                        raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
+                forcing_data = pd.read_csv(os.path.join(forcing_dir, var_lst[k] + ".csv"))
+                df_date = forcing_data[['year', 'month', 'day']]
+                date = pd.to_datetime(df_date).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                chosen_data = forcing_data[gage_id_lst].values[ind1, :]
+                x[:, ind2, k] = chosen_data.T
+        elif self.region == "BR":
+            for j in range(len(var_lst)):
+                for k in range(len(gage_id_lst)):
+                    data_obs = self.read_br_basin_forcing(gage_id_lst[k], t_range, var_lst[j])
+                    x[k, :, j] = data_obs
+        elif self.region == "CL":
+            for k in range(len(var_lst)):
+                for tmp in os.listdir(self.data_source_description["CAMELS_DIR"]):
+                    if fnmatch.fnmatch(tmp, '*' + var_lst[k]):
+                        tmp_ = os.path.join(self.data_source_description["CAMELS_DIR"], tmp)
+                        if os.path.isdir(tmp_):
+                            forcing_file = os.path.join(tmp_, os.listdir(tmp_)[0])
+                forcing_data = pd.read_csv(forcing_file, sep="\t", index_col=0)
+                date = pd.to_datetime(forcing_data.index.values).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                station_ids = self.read_object_ids()
+                assert (all(x < y for x, y in zip(station_ids, station_ids[1:])))
+                [s, ind3, ind4] = np.intersect1d(station_ids, gage_id_lst, return_indices=True)
+                chosen_data = forcing_data.iloc[ind1, ind3].replace("\s+", np.nan, regex=True)
+                x[:, ind2, k] = chosen_data.values.T
+        elif self.region == "GB":
+            for j in range(len(var_lst)):
+                for k in range(len(gage_id_lst)):
+                    data_forcing = self.read_gb_gage_flow_forcing(gage_id_lst[k], t_range, var_lst[j])
+                    x[k, :, j] = data_forcing
+        elif self.region == "YR":
+            for k in range(len(gage_id_lst)):
+                forcing_file = os.path.join(self.data_source_description["CAMELS_FORCING_DIR"], gage_id_lst[k],
+                                            "forcing.csv")
+                forcing_data = pd.read_csv(forcing_file, sep=",")
+                date = pd.to_datetime(forcing_data["date"]).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                for j in range(len(var_lst)):
+                    x[k, ind2, j] = forcing_data[var_lst[j]].values[ind1]
+        elif self.region == "CA":
+            for k in range(len(gage_id_lst)):
+                canopex_id = self.camels_sites[
+                    self.camels_sites["STATION_ID"] == "'" + gage_id_lst[k] + "'"]["CANOPEX_ID"].values[0]
+                forcing_file = os.path.join(self.data_source_description["CAMELS_FLOW_DIR"], str(canopex_id) + ".dly")
+                read_forcing_file = pd.read_csv(forcing_file, header=None).values.tolist()
+
+                forcing_date = []
+                for j in range(len(var_lst)):
+                    forcing_data = []
+                    for one_site in read_forcing_file:
+                        forcing_date.append(hydro_utils.t2dt(int(one_site[0][:8].replace(" ", "0"))))
+                        all_data = one_site[0].split(" ")
+                        real_data = [one_data for one_data in all_data if one_data != ""]
+                        if var_lst[j] == "prcp":
+                            forcing_data.append(float(real_data[-5]))
+                        elif var_lst[j] == "tmax":
+                            forcing_data.append(float(real_data[-2]))
+                        elif var_lst[j] == "tmin":
+                            forcing_data.append(float(real_data[-1]))
+                        else:
+                            raise NotImplementedError("No such forcing type in CANOPEX now!")
+                    date = pd.to_datetime(forcing_date).values.astype('datetime64[D]')
+                    [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                    x[k, ind2, j] = np.array(forcing_data)[ind1]
+        elif self.region == "CE":
+            for k in range(len(gage_id_lst)):
+                forcing_file = os.path.join(self.data_source_description["CAMELS_FORCING_DIR"],
+                                            "ID_" + str(gage_id_lst[k]) + ".csv")
+                forcing_data = pd.read_csv(forcing_file, sep=";")
+                df_date = forcing_data[["YYYY", "MM", "DD"]]
+                df_date.columns = ['year', 'month', 'day']
+                date = pd.to_datetime(df_date).values.astype('datetime64[D]')
+                [c, ind1, ind2] = np.intersect1d(date, t_range_list, return_indices=True)
+                for j in range(len(var_lst)):
+                    x[k, ind2, j] = forcing_data[var_lst[j]].values[ind1]
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
         return x
 
     def read_attr_all(self):
         data_folder = self.data_source_description["CAMELS_ATTR_DIR"]
+        key_lst = self.data_source_description["CAMELS_ATTR_KEY_LST"]
         f_dict = dict()  # factorize dict
         var_dict = dict()
         var_lst = list()
         out_lst = list()
-        key_lst = ['topo', 'clim', 'hydro', 'vege', 'soil', 'geol']
         gage_dict = self.camels_sites
+        if self.region == "US":
+            camels_str = "camels_"
+            sep_ = ';'
+        elif self.region == "BR":
+            camels_str = "camels_br_"
+            sep_ = '\s+'
+        elif self.region == "GB":
+            camels_str = "CAMELS_GB_"
+            sep_ = ','
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
         for key in key_lst:
-            data_file = os.path.join(data_folder, 'camels_' + key + '.txt')
-            data_temp = pd.read_csv(data_file, sep=';')
+            data_file = os.path.join(data_folder, camels_str + key + '.txt')
+            if self.region == "GB":
+                data_file = os.path.join(data_folder, camels_str + key + '_attributes.csv')
+            data_temp = pd.read_csv(data_file, sep=sep_)
             var_lst_temp = list(data_temp.columns[1:])
             var_dict[key] = var_lst_temp
             var_lst.extend(var_lst_temp)
@@ -190,14 +976,116 @@ def read_attr_all(self):
         out = np.concatenate(out_lst, 1)
         return out, var_lst, var_dict, f_dict
 
-    def read_constant_cols(self, usgs_id_lst=None, var_lst=None, is_return_dict=False):
-        attr_all, var_lst_all, var_dict, f_dict = self.read_attr_all()
+    def read_attr_all_in_one_file(self):
+        """
+        Read all attr data in CAMELS_AUS or CAMELS_CL
+
+        Returns
+        -------
+        np.array
+            all attr data in CAMELS_AUS or CAMELS_CL
+        """
+        if self.region == "AUS":
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_DIR"],
+                                         "CAMELS_AUS_Attributes-Indices_MasterTable.csv")
+            all_attr = pd.read_csv(attr_all_file, sep=',')
+        elif self.region == "CL":
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_ATTR_DIR"], "1_CAMELScl_attributes.txt")
+            all_attr_tmp = pd.read_csv(attr_all_file, sep='\t', index_col=0)
+            all_attr = pd.DataFrame(all_attr_tmp.values.T, index=all_attr_tmp.columns, columns=all_attr_tmp.index)
+        elif self.region == "CA":
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_ATTR_DIR"],
+                                         "HYSETS_watershed_properties.txt")
+            all_attr_tmp = pd.read_csv(attr_all_file, sep=';', index_col=0)
+            all_attr = all_attr_tmp[all_attr_tmp["Official_ID"].isin(self.read_object_ids())]
+        elif self.region == "CE":
+            attr_all_file = os.path.join(self.data_source_description["CAMELS_ATTR_DIR"], "Catchment_attributes.csv")
+            all_attr = pd.read_csv(attr_all_file, sep=';')
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
+        # gage_all_attr = all_attr[all_attr['station_id'].isin(gage_id_lst)]
+        var_lst = self.get_constant_cols().tolist()
+        data_temp = all_attr[var_lst]
+        # for factorized data, we need factorize all gages' data to keep the factorized number same all the time
+        n_gage = len(self.read_object_ids())
+        out_temp = np.full([n_gage, len(var_lst)], np.nan)
+        f_dict = {}
+        out_lst = []
+        k = 0
+        for field in var_lst:
+            if is_string_dtype(data_temp[field]):
+                value, ref = pd.factorize(data_temp[field], sort=True)
+                out_temp[:, k] = value
+                f_dict[field] = ref.tolist()
+            elif is_numeric_dtype(data_temp[field]):
+                out_temp[:, k] = data_temp[field].values
+            k = k + 1
+            out_lst.append(out_temp)
+        out = np.concatenate(out_lst, 1)
+        # keep same format with CAMELS_US
+        return out, var_lst, None, f_dict
+
+    def read_attr_all_yr(self):
+        var_lst = self.get_constant_cols().tolist()
+        gage_id_lst = self.read_object_ids()
+        # for factorized data, we need factorize all gages' data to keep the factorized number same all the time
+        n_gage = len(self.read_object_ids())
+        c = np.empty([n_gage, len(var_lst)], dtype=object)
+        for k in range(n_gage):
+            attr_file = os.path.join(self.data_source_description["CAMELS_ATTR_DIR"], gage_id_lst[k],
+                                     "attributes.json")
+            attr_data = hydro_utils.unserialize_json_ordered(attr_file)
+            for j in range(len(var_lst)):
+                c[k, j] = attr_data[var_lst[j]]
+        data_temp = pd.DataFrame(c, columns=var_lst)
+        out_temp = np.full([n_gage, len(var_lst)], np.nan)
+        f_dict = {}
+        k = 0
+        for field in var_lst:
+            if field in ["high_prec_timing", "low_prec_timing"]:
+                # string type
+                value, ref = pd.factorize(data_temp[field], sort=True)
+                out_temp[:, k] = value
+                f_dict[field] = ref.tolist()
+            else:
+                out_temp[:, k] = data_temp[field].values
+            k = k + 1
+        # keep same format with CAMELS_US
+        return out_temp, var_lst, None, f_dict
+
+    def read_constant_cols(self, gage_id_lst=None, var_lst=None, is_return_dict=False) -> Union[tuple, np.array]:
+        """
+        Read Attributes data
+
+        Parameters
+        ----------
+        gage_id_lst
+            station ids
+        var_lst
+            attribute variable types
+        is_return_dict
+            if true, return var_dict and f_dict for CAMELS_US
+        Returns
+        -------
+        Union[tuple, np.array]
+            if attr var type is str, return factorized data.
+            When we need to know what a factorized value represents, we need return a tuple;
+            otherwise just return an array
+        """
+        if self.region in ["BR", "GB", "US"]:
+            attr_all, var_lst_all, var_dict, f_dict = self.read_attr_all()
+        elif self.region in ["AUS", "CA", "CE", "CL"]:
+            attr_all, var_lst_all, var_dict, f_dict = self.read_attr_all_in_one_file()
+        elif self.region == "YR":
+            attr_all, var_lst_all, var_dict, f_dict = self.read_attr_all_yr()
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
         ind_var = list()
         for var in var_lst:
             ind_var.append(var_lst_all.index(var))
-        gage_dict = self.camels_sites
-        id_lst_all = gage_dict['gauge_id'].values
-        c, ind_grid, ind2 = np.intersect1d(id_lst_all, usgs_id_lst, return_indices=True)
+        id_lst_all = self.read_object_ids()
+        # TODO: Notice the sequence of station ids !!!!!!
+        c, ind_grid, ind2 = np.intersect1d(id_lst_all, gage_id_lst, return_indices=True)
         temp = attr_all[ind_grid, :]
         out = temp[:, ind_var]
         if is_return_dict:
@@ -205,6 +1093,46 @@ def read_constant_cols(self, usgs_id_lst=None, var_lst=None, is_return_dict=Fals
         else:
             return out
 
+    def read_basin_area(self, object_ids) -> np.array:
+        if self.region == "US":
+            return self.read_constant_cols(object_ids, ['area_gages2'], is_return_dict=False)
+        elif self.region == "AUS":
+            return self.read_constant_cols(object_ids, ['catchment_area'], is_return_dict=False)
+        elif self.region in ["BR", "CL", "GB", "YR"]:
+            return self.read_constant_cols(object_ids, ['area'], is_return_dict=False)
+        elif self.region == "CA":
+            return self.read_constant_cols(object_ids, ['Drainage_Area_km2'], is_return_dict=False)
+        elif self.region == "CE":
+            return self.read_constant_cols(object_ids, ['area_calc'], is_return_dict=False)
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
+
+    def read_mean_prep(self, object_ids) -> np.array:
+        if self.region in ["US", "AUS", "BR", "GB", "CE"]:
+            return self.read_constant_cols(object_ids, ['p_mean'], is_return_dict=False)
+        elif self.region == "CL":
+            # there are different p_mean values for different forcings, here we chose p_mean_cr2met now
+            return self.read_constant_cols(object_ids, ['p_mean_cr2met'], is_return_dict=False)
+        elif self.region == "YR":
+            return self.read_constant_cols(object_ids, ['pre_mean'], is_return_dict=False)
+        elif self.region == "CA":
+            # There is no p_mean attr, hence we have to calculate from forcing data directly
+            prcp_means = []
+            for k in range(len(object_ids)):
+                canopex_id = self.camels_sites[
+                    self.camels_sites["STATION_ID"] == "'" + object_ids[k] + "'"]["CANOPEX_ID"].values[0]
+                forcing_file = os.path.join(self.data_source_description["CAMELS_FLOW_DIR"], str(canopex_id) + ".dly")
+                read_forcing_file = pd.read_csv(forcing_file, header=None).values.tolist()
+                prcp_data = []
+                for one_site in read_forcing_file:
+                    all_data = one_site[0].split(" ")
+                    real_data = [one_data for one_data in all_data if one_data != ""]
+                    prcp_data.append(float(real_data[-5]))
+                prcp_means.append(np.mean(np.array(prcp_data)))
+            return np.array(prcp_means)
+        else:
+            raise NotImplementedError(CAMELS_NO_DATASET_ERROR_LOG)
+
     def read_other_cols(self, object_ids=None, other_cols: dict = None, **kwargs):
         # TODO: FDC for test period should keep same with that in training period
         out_dict = {}


Test patch:
diff --git a/test/test_camels.py b/test/test_camels.py
index 9af2b81..e6b40d3 100644
--- a/test/test_camels.py
+++ b/test/test_camels.py
@@ -1,20 +1,277 @@
 import os
 import unittest
+
+import numpy as np
+
 import definitions
 from hydrobench.data.data_camels import Camels
 
 
 class CamelsTestCase(unittest.TestCase):
     def setUp(self) -> None:
-        self.camels_path = os.path.join(definitions.DATASET_DIR, "camels")
+        self.camels_aus_path = os.path.join(definitions.DATASET_DIR, "camels", "camels_aus")
+        self.camels_br_path = os.path.join(definitions.DATASET_DIR, "camels", "camels_br")
+        self.camels_cl_path = os.path.join(definitions.DATASET_DIR, "camels", "camels_cl")
+        self.camels_gb_path = os.path.join(definitions.DATASET_DIR, "camels", "camels_gb")
+        self.camels_us_path = os.path.join(definitions.DATASET_DIR, "camels", "camels_us")
+        self.camels_yr_path = os.path.join(definitions.DATASET_DIR, "camels", "camels_yr")
+        self.canopex_path = os.path.join(definitions.DATASET_DIR, "canopex")
+        self.lamah_ce_path = os.path.join(definitions.DATASET_DIR, "lamah_ce")
+        self.aus_region = "AUS"
+        self.br_region = "BR"
+        self.cl_region = "CL"
+        self.gb_region = "GB"
+        self.us_region = "US"
+        self.yr_region = "YR"
+        self.ca_region = "CA"
+        self.lamah_ce_region = "CE"
 
     def test_download_camels(self):
-        camels = Camels(self.camels_path, download=True)
-        self.assertTrue(os.path.isfile(os.path.join(self.camels_path, "basin_set_full_res", "HCDN_nhru_final_671.shp")))
+        camels_us = Camels(self.camels_us_path, download=True)
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.camels_us_path, "basin_set_full_res", "HCDN_nhru_final_671.shp")))
+
+    def test_read_camels_us(self):
+        camels_us = Camels(self.camels_us_path, download=False, region=self.us_region)
+        gage_ids = camels_us.read_object_ids()
+        self.assertEqual(gage_ids.size, 671)
+        attrs = camels_us.read_constant_cols(gage_ids[:5], var_lst=["soil_conductivity", "elev_mean", "geol_1st_class"])
+        np.testing.assert_almost_equal(attrs, np.array(
+            [[1.10652248, 250.31, 10.], [2.37500506, 92.68, 0.], [1.28980735, 143.8, 10.],
+             [1.37329168, 247.8, 10.], [2.61515428, 310.38, 7.]]))
+        forcings = camels_us.read_relevant_cols(gage_ids[:5], ["1990-01-01", "2010-01-01"],
+                                                var_lst=["dayl", "prcp", "srad"])
+        np.testing.assert_array_equal(forcings.shape, np.array([5, 7305, 3]))
+        flows = camels_us.read_target_cols(gage_ids[:5], ["1990-01-01", "2010-01-01"], target_cols=["usgsFlow"])
+        np.testing.assert_array_equal(flows.shape, np.array([5, 7305, 1]))
+        streamflow_types = camels_us.get_target_cols()
+        np.testing.assert_array_equal(streamflow_types, np.array(["usgsFlow"]))
+        focing_types = camels_us.get_relevant_cols()
+        np.testing.assert_array_equal(focing_types, np.array(['dayl', 'prcp', 'srad', 'swe', 'tmax', 'tmin', 'vp']))
+        attr_types = camels_us.get_constant_cols()
+        np.testing.assert_array_equal(attr_types[:3], np.array(['gauge_lat', 'gauge_lon', 'elev_mean']))
+
+    def test_download_camels_aus(self):
+        camels_aus = Camels(self.camels_aus_path, download=True, region="AUS")
+        self.assertTrue(os.path.isfile(os.path.join(self.camels_aus_path, "05_hydrometeorology", "05_hydrometeorology",
+                                                    "01_precipitation_timeseries", "precipitation_AWAP.csv")))
+
+    def test_read_camels_aus(self):
+        camels_aus = Camels(self.camels_aus_path, download=False, region=self.aus_region)
+        gage_ids = camels_aus.read_object_ids()
+        self.assertEqual(gage_ids.size, 222)
+        attrs = camels_aus.read_constant_cols(gage_ids[:5], var_lst=["catchment_area", "slope_fdc", "geol_sec"])
+        np.testing.assert_array_equal(attrs, np.array(
+            [[1.25773e+04, 3.66793e-01, 6.00000e+00], [1.13929e+04, 3.29998e-01, 6.00000e+00],
+             [5.65300e+02, 1.78540e-02, 6.00000e+00], [4.58300e+02, 5.00234e-01, 6.00000e+00],
+             [7.73170e+03, 3.74751e-01, 1.00000e+00]]))
+        forcings = camels_aus.read_relevant_cols(gage_ids[:5], ["1990-01-01", "2010-01-01"],
+                                                 var_lst=["precipitation_AWAP", "et_morton_actual_SILO", "tmin_SILO"])
+        np.testing.assert_array_equal(forcings.shape, np.array([5, 7305, 3]))
+        flows = camels_aus.read_target_cols(gage_ids[:5], ["1990-01-01", "2010-01-01"],
+                                            target_cols=["streamflow_MLd", "streamflow_mmd"])
+        np.testing.assert_array_equal(flows.shape, np.array([5, 7305, 2]))
+        streamflow_types = camels_aus.get_target_cols()
+        np.testing.assert_array_equal(streamflow_types, np.array(
+            ["streamflow_MLd", "streamflow_MLd_inclInfilled", "streamflow_mmd", "streamflow_QualityCodes"]))
+        focing_types = camels_aus.get_relevant_cols()
+        np.testing.assert_array_equal(focing_types, np.array(
+            ['precipitation_AWAP', 'precipitation_SILO', 'precipitation_var_AWAP', 'et_morton_actual_SILO',
+             'et_morton_point_SILO', 'et_morton_wet_SILO', 'et_short_crop_SILO', 'et_tall_crop_SILO',
+             'evap_morton_lake_SILO', 'evap_pan_SILO', 'evap_syn_SILO', 'solarrad_AWAP', 'tmax_AWAP', 'tmin_AWAP',
+             'vprp_AWAP', 'mslp_SILO', 'radiation_SILO', 'rh_tmax_SILO', 'rh_tmin_SILO', 'tmax_SILO', 'tmin_SILO',
+             'vp_deficit_SILO', 'vp_SILO']))
+        attr_types = camels_aus.get_constant_cols()
+        np.testing.assert_array_equal(attr_types[:3],
+                                      np.array(['station_name', 'drainage_division', 'river_region']))
+
+    def test_download_camels_br(self):
+        camels_br = Camels(self.camels_br_path, download=True, region="BR")
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.camels_br_path, "01_CAMELS_BR_attributes", "01_CAMELS_BR_attributes",
+                                        "CAMELS_BR_attributes_description.xlsx")))
+
+    def test_read_camels_br(self):
+        camels_br = Camels(self.camels_br_path, download=False, region=self.br_region)
+        gage_ids = camels_br.read_object_ids()
+        self.assertEqual(gage_ids.size, 897)
+        attrs = camels_br.read_constant_cols(gage_ids[:5], var_lst=["geol_class_1st", "p_mean", "runoff_ratio"])
+        np.testing.assert_array_equal(attrs, np.array(
+            [[8., 6.51179, 0.55336], [6., 5.38941, 0.72594], [6., 5.70191, 0.759], [6., 5.19877, 0.39463],
+             [8., 5.49805, 0.38579]]))
+        forcings = camels_br.read_relevant_cols(gage_ids[:5], ["1995-01-01", "2015-01-01"],
+                                                var_lst=["precipitation_chirps", "evapotransp_gleam",
+                                                         "potential_evapotransp_gleam", "temperature_min_cpc"])
+        np.testing.assert_array_equal(forcings.shape, np.array([5, 7305, 4]))
+        # start from 1995/01/01 end with 2017/04/30
+        flows = camels_br.read_target_cols(gage_ids[:5], ["1995-01-01", "2015-01-01"],
+                                           target_cols=["streamflow_m3s", "streamflow_mm_selected_catchments"])
+        np.testing.assert_array_equal(flows.shape, np.array([5, 7305, 2]))
+        streamflow_types = camels_br.get_target_cols()
+        np.testing.assert_array_equal(streamflow_types, np.array(
+            ['streamflow_m3s', 'streamflow_mm_selected_catchments', 'streamflow_simulated']))
+        focing_types = camels_br.get_relevant_cols()
+        np.testing.assert_array_equal(focing_types, np.array(
+            ['precipitation_chirps', 'precipitation_mswep', 'precipitation_cpc', 'evapotransp_gleam', 'evapotransp_mgb',
+             'potential_evapotransp_gleam', 'temperature_min_cpc', 'temperature_mean_cpc', 'temperature_max_cpc']))
+        attr_types = camels_br.get_constant_cols()
+        np.testing.assert_array_equal(attr_types[:3], np.array(['p_mean', 'pet_mean', 'et_mean']))
+
+    def test_download_camels_cl(self):
+        camels_cl = Camels(self.camels_cl_path, download=True, region="CL")
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.camels_cl_path, "1_CAMELScl_attributes", "1_CAMELScl_attributes.txt")))
+
+    def test_read_camels_cl(self):
+        camels_cl = Camels(self.camels_cl_path, download=False, region=self.cl_region)
+        gage_ids = camels_cl.read_object_ids()
+        self.assertEqual(gage_ids.size, 516)
+        # TODO: NaN value
+        attrs = camels_cl.read_constant_cols(gage_ids[:5], var_lst=["geol_class_1st", "crop_frac"])
+        np.testing.assert_array_equal(attrs, np.array([[9., 0.], [9., 53.], [9., 64.], [10., 144.], [10., 104.]]))
+        # TODO: NaN value, mainly for pet_8d_modis
+        forcings = camels_cl.read_relevant_cols(gage_ids[:5], ["1995-01-01", "2015-01-01"],
+                                                var_lst=["pet_8d_modis", "precip_cr2met", "swe"])
+        np.testing.assert_array_equal(forcings.shape, np.array([5, 7305, 3]))
+        flows = camels_cl.read_target_cols(gage_ids[:5], ["1995-01-01", "2015-01-01"],
+                                           target_cols=["streamflow_m3s", "streamflow_mm"])
+        np.testing.assert_array_equal(flows.shape, np.array([5, 7305, 2]))
+        streamflow_types = camels_cl.get_target_cols()
+        np.testing.assert_array_equal(streamflow_types, np.array(['streamflow_m3s', 'streamflow_mm']))
+        focing_types = camels_cl.get_relevant_cols()
+        np.testing.assert_array_equal(focing_types, np.array(
+            ['precip_cr2met', 'precip_chirps', 'precip_mswep', 'precip_tmpa', 'tmin_cr2met', 'tmax_cr2met',
+             'tmean_cr2met', 'pet_8d_modis', 'pet_hargreaves', 'swe']))
+        attr_types = camels_cl.get_constant_cols()
+        np.testing.assert_array_equal(attr_types[:3], np.array(['gauge_name', 'gauge_lat', 'gauge_lon']))
+
+    def test_download_camels_gb(self):
+        camels_gb = Camels(self.camels_gb_path, download=True, region=self.gb_region)
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.camels_gb_path, "8344e4f3-d2ea-44f5-8afa-86d2987543a9",
+                                        "8344e4f3-d2ea-44f5-8afa-86d2987543a9", "data",
+                                        "CAMELS_GB_climatic_attributes.csv")))
+
+    def test_read_camels_gb(self):
+        camels_gb = Camels(self.camels_gb_path, download=False, region=self.gb_region)
+        gage_ids = camels_gb.read_object_ids()
+        self.assertEqual(gage_ids.size, 671)
+        attrs = camels_gb.read_constant_cols(gage_ids[:5], var_lst=["p_mean", "slope_fdc", "gauge_name"])
+        np.testing.assert_array_equal(attrs, np.array(
+            [[2.29, 1.94, 596.], [2.31, 1.95, 670.], [2.65, 4.01, 647.], [2.31, 1.54, 393.], [2.29, 1.47, 217.]]))
+        forcings = camels_gb.read_relevant_cols(gage_ids[:5], ["1995-01-01", "2015-01-01"],
+                                                var_lst=["precipitation", "pet", "temperature", "peti"])
+        np.testing.assert_array_equal(forcings.shape, np.array([5, 7305, 4]))
+        flows = camels_gb.read_target_cols(gage_ids[:5], ["1995-01-01", "2015-01-01"],
+                                           target_cols=["discharge_spec", "discharge_vol"])
+        np.testing.assert_array_equal(flows.shape, np.array([5, 7305, 2]))
+        streamflow_types = camels_gb.get_target_cols()
+        np.testing.assert_array_equal(streamflow_types, np.array(["discharge_spec", "discharge_vol"]))
+        focing_types = camels_gb.get_relevant_cols()
+        np.testing.assert_array_equal(focing_types, np.array(
+            ["precipitation", "pet", "temperature", "peti", "humidity", "shortwave_rad", "longwave_rad", "windspeed"]))
+        attr_types = camels_gb.get_constant_cols()
+        np.testing.assert_array_equal(attr_types[:3], np.array(["p_mean", "pet_mean", "aridity"]))
+
+    def test_download_camels_yr(self):
+        camels_yr = Camels(self.camels_yr_path, download=True, region=self.yr_region)
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.camels_yr_path, "8344e4f3-d2ea-44f5-8afa-86d2987543a9",
+                                        "8344e4f3-d2ea-44f5-8afa-86d2987543a9", "data",
+                                        "CAMELS_GB_climatic_attributes.csv")))
+
+    def test_read_camels_yr(self):
+        camels_yr = Camels(self.camels_yr_path, download=False, region=self.yr_region)
+        gage_ids = camels_yr.read_object_ids()
+        self.assertEqual(gage_ids.size, 102)
+        attrs = camels_yr.read_constant_cols(gage_ids[:5], var_lst=["area", "barren", "bdticm"])
+        np.testing.assert_almost_equal(attrs, np.array(
+            [[3.11520000e+04, 2.98706264e-03, 3.33904449e+03], [4.27056000e+05, 2.30162622e-02, 1.80570119e+03],
+             [3.80128000e+05, 2.47549979e-02, 1.77874628e+03], [7.33561000e+05, 1.41340180e-02, 2.01143843e+03],
+             [2.04213000e+05, 7.75394506e-03, 1.53321208e+03]]), decimal=3)
+        forcings = camels_yr.read_relevant_cols(gage_ids[:5], ["1995-01-01", "2015-01-01"],
+                                                var_lst=["pre", "evp", "gst_mean", "prs_mean"])
+        np.testing.assert_array_equal(forcings.shape, np.array([5, 7305, 4]))
+        flows = camels_yr.read_target_cols(gage_ids[:5], ["1995-01-01", "2015-01-01"], target_cols=["normalized_q"])
+        np.testing.assert_array_equal(flows.shape, np.array([5, 7305, 1]))
+        streamflow_types = camels_yr.get_target_cols()
+        np.testing.assert_array_equal(streamflow_types, np.array(["normalized_q"]))
+        focing_types = camels_yr.get_relevant_cols()
+        np.testing.assert_array_equal(focing_types, np.array(
+            ["pre", "evp", "gst_mean", "prs_mean", "tem_mean", "rhu", "win_mean", "gst_min", "prs_min", "tem_min",
+             "gst_max", "prs_max", "tem_max", "ssd", "win_max"]))
+        attr_types = camels_yr.get_constant_cols()
+        np.testing.assert_array_equal(attr_types[:3], np.array(["area", "barren", "bdticm"]))
+
+    def test_download_canopex(self):
+        canopex = Camels(self.canopex_path, download=True, region=self.ca_region)
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.canopex_path, "CANOPEX_NRCAN_ASCII", "CANOPEX_NRCAN_ASCII", "1.dly")))
+
+    def test_read_canopex_data(self):
+        canopex = Camels(self.canopex_path, download=False, region=self.ca_region)
+        gage_ids = canopex.read_object_ids()
+        self.assertEqual(gage_ids.size, 611)
+        attrs = canopex.read_constant_cols(gage_ids[:5],
+                                           var_lst=["Drainage_Area_km2", "Land_Use_Grass_frac", "Permeability_logk_m2"])
+        np.testing.assert_almost_equal(attrs, np.array(
+            [[3.28438700e+02, 6.94000000e-02, - 1.35251586e+01], [3.43515600e+02, 6.16000000e-02, - 1.42367348e+01],
+             [1.45554950e+03, 3.59000000e-02, - 1.50071080e+01], [5.64820300e+02, 3.05000000e-02, - 1.51002546e+01],
+             [1.05383090e+03, 3.27000000e-02, - 1.51999998e+01]]), decimal=3)
+        forcings = canopex.read_relevant_cols(gage_ids[:5], ["1990-01-01", "2010-01-01"],
+                                              var_lst=["prcp", "tmax", "tmin"])
+        np.testing.assert_array_equal(forcings.shape, np.array([5, 7305, 3]))
+        flows = canopex.read_target_cols(gage_ids[:5], ["1990-01-01", "2010-01-01"], target_cols=["discharge"])
+        np.testing.assert_array_equal(flows.shape, np.array([5, 7305, 1]))
+        streamflow_types = canopex.get_target_cols()
+        np.testing.assert_array_equal(streamflow_types, np.array(["discharge"]))
+        focing_types = canopex.get_relevant_cols()
+        np.testing.assert_array_equal(focing_types, np.array(["prcp", "tmax", "tmin"]))
+        attr_types = canopex.get_constant_cols()
+        np.testing.assert_array_equal(attr_types[:3], np.array(["Source", "Name", "Official_ID"]))
+
+    def test_download_lamah_ce(self):
+        lamah_ce = Camels(self.lamah_ce_path, download=True, region=self.lamah_ce_region)
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.lamah_ce_path, "CANOPEX_NRCAN_ASCII", "CANOPEX_NRCAN_ASCII", "1.dly")))
+
+    def test_read_lamah_ce(self):
+        lamah_ce = Camels(self.lamah_ce_path, download=False, region=self.lamah_ce_region)
+        gage_ids = lamah_ce.read_object_ids()
+        self.assertEqual(gage_ids.size, 859)
+        attrs = lamah_ce.read_constant_cols(gage_ids[:5], var_lst=['area_calc', 'elev_mean', 'elev_med'])
+        np.testing.assert_almost_equal(attrs, np.array(
+            [[4668.379, 1875., 1963.], [102.287, 1775., 1827.], [536.299, 1844., 1916.], [66.286, 1894., 1907.],
+             [72.448, 1774., 1796.]]), decimal=3)
+        forcings = lamah_ce.read_relevant_cols(gage_ids[:5], ["1990-01-01", "2010-01-01"],
+                                               var_lst=["2m_temp_max", "prec", "volsw_4"])
+        np.testing.assert_array_equal(forcings.shape, np.array([5, 7305, 3]))
+        flows = lamah_ce.read_target_cols(gage_ids[:5], ["1990-01-01", "2010-01-01"], target_cols=["qobs"])
+        np.testing.assert_array_equal(flows.shape, np.array([5, 7305, 1]))
+        streamflow_types = lamah_ce.get_target_cols()
+        np.testing.assert_array_equal(streamflow_types, np.array(["qobs"]))
+        focing_types = lamah_ce.get_relevant_cols()
+        np.testing.assert_array_equal(focing_types, np.array(
+            ["2m_temp_max", "2m_temp_mean", "2m_temp_min", "2m_dp_temp_max", "2m_dp_temp_mean",
+             "2m_dp_temp_min", "10m_wind_u", "10m_wind_v", "fcst_alb", "lai_high_veg", "lai_low_veg",
+             "swe", "surf_net_solar_rad_max", "surf_net_solar_rad_mean", "surf_net_therm_rad_max",
+             "surf_net_therm_rad_mean", "surf_press", "total_et", "prec", "volsw_123", "volsw_4"]))
+        attr_types = lamah_ce.get_constant_cols()
+        np.testing.assert_array_equal(attr_types, np.array(
+            ['area_calc', 'elev_mean', 'elev_med', 'elev_std', 'elev_ran', 'slope_mean', 'mvert_dist', 'mvert_ang',
+             'elon_ratio', 'strm_dens', 'p_mean', 'et0_mean', 'eta_mean', 'arid_1', 'arid_2', 'p_season', 'frac_snow',
+             'hi_prec_fr', 'hi_prec_du', 'hi_prec_ti', 'lo_prec_fr', 'lo_prec_du', 'lo_prec_ti', 'lc_dom', 'agr_fra',
+             'bare_fra', 'forest_fra', 'glac_fra', 'lake_fra', 'urban_fra', 'lai_max', 'lai_diff', 'ndvi_max',
+             'ndvi_min', 'gvf_max', 'gvf_diff', 'bedrk_dep', 'root_dep', 'soil_poros', 'soil_condu', 'soil_tawc',
+             'sand_fra', 'silt_fra', 'clay_fra', 'grav_fra', 'oc_fra', 'gc_dom', 'gc_ig_fra', 'gc_mt_fra', 'gc_pa_fra',
+             'gc_pb_fra', 'gc_pi_fra', 'gc_py_fra', 'gc_sc_fra', 'gc_sm_fra', 'gc_ss_fra', 'gc_su_fra', 'gc_va_fra',
+             'gc_vb_fra', 'gc_wb_fra', 'geol_perme', 'geol_poros']))
 
-    def test_download_camels_after_downloading(self):
-        camels = Camels(self.camels_path, download=True)
-        self.assertTrue(os.path.isfile(os.path.join(self.camels_path, "basin_set_full_res", "HCDN_nhru_final_671.shp")))
+    def test_ca_p_mean(self):
+        canopex = Camels(self.canopex_path, download=False, region=self.ca_region)
+        gage_ids = canopex.read_object_ids()
+        p_mean = canopex.read_mean_prep(gage_ids[:5])
+        np.testing.assert_almost_equal(p_mean, np.array([2.91712073, 3.14145935, 3.12958083, 3.09248435, 3.04431583]))
 
 
 if __name__ == '__main__':