From dacdedec035bf8b0f663983d09cbec9ce8654151 Mon Sep 17 00:00:00 2001
From: annajungbluth <anna.jungbluth1993@gmail.com>
Date: Fri, 15 Mar 2024 11:06:01 +0000
Subject: [PATCH] added todos, deleted duplicate functions

---
 notebooks/goes-scripts/goes_downsample.py    |  24 -
 notebooks/goes-scripts/goes_l1b_download.py  | 649 -------------------
 notebooks/goes-scripts/xy_to_latlon.py       |  73 ---
 rs_tools/_src/geoprocessing/grid.py          |   2 +-
 rs_tools/_src/geoprocessing/modis/interp.py  |   1 +
 rs_tools/_src/geoprocessing/modis/rescale.py |   1 +
 rs_tools/_src/geoprocessing/reproject.py     |   1 +
 rs_tools/_src/modules/__init__.py            |   0
 rs_tools/_src/modules/modis/__init__.py      |   0
 rs_tools/_src/modules/modis/download.py      |   0
 rs_tools/_src/preprocessing/normalize.py     |   6 +-
 rs_tools/_src/utils/crs_ops.py               | 120 ----
 rs_tools/_src/utils/math.py                  |   1 +
 13 files changed, 8 insertions(+), 870 deletions(-)
 delete mode 100644 notebooks/goes-scripts/goes_downsample.py
 delete mode 100644 notebooks/goes-scripts/goes_l1b_download.py
 delete mode 100644 notebooks/goes-scripts/xy_to_latlon.py
 delete mode 100644 rs_tools/_src/modules/__init__.py
 delete mode 100644 rs_tools/_src/modules/modis/__init__.py
 delete mode 100644 rs_tools/_src/modules/modis/download.py
 delete mode 100644 rs_tools/_src/utils/crs_ops.py

diff --git a/notebooks/goes-scripts/goes_downsample.py b/notebooks/goes-scripts/goes_downsample.py
deleted file mode 100644
index 41c3ab6..0000000
--- a/notebooks/goes-scripts/goes_downsample.py
+++ /dev/null
@@ -1,24 +0,0 @@
-from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator
-import numpy as np
-from xy_to_latlon import calc_latlon
-
-def downsample(ds, target_points):
-    """
-    target_points
-    """
-    # attach latitude & longitude coordinates to xr.Datasets
-    ds = calc_latlon(ds)
-    
-    # extract 1d arrays of latitudes and longitudes
-    lat = ds.lat.to_numpy().flatten()
-    lon = ds.lon.to_numpy().flatten()
-    
-    # turn into 2d array of latitudes and longitudes
-    points = np.vstack((lon, lat)).T
-    
-    # initialise interpolation
-    nn_interpolation = NearestNDInterpolator(points, ds.Rad)
-    interpolated_nn = nn_interpolation(target_points)
-    
-    # create new xr.Dataset with lowres data 
-    ...
\ No newline at end of file
diff --git a/notebooks/goes-scripts/goes_l1b_download.py b/notebooks/goes-scripts/goes_l1b_download.py
deleted file mode 100644
index 48b8980..0000000
--- a/notebooks/goes-scripts/goes_l1b_download.py
+++ /dev/null
@@ -1,649 +0,0 @@
-# Builtin modules
-import os
-import shutil
-import subprocess
-import warnings
-from datetime import datetime, timedelta
-
-# External modules
-from google.cloud import storage
-from dateutil.parser import parse as parse_date
-import numpy as np
-import xarray as xr
-
-# set path to google cloud credentials
-os.environ["GOOGLE_APPLICATION_CREDENTIALS"] = 'ADD GOOGLE CREDENTIALS FILE HERE'
-
-
-if "GOOGLE_APPLICATION_CREDENTIALS" not in os.environ:
-    raise ImportError(
-        """
-        'GOOGLE_APPLICATION_CREDENTIALS' is not set, this is required for IO
-        operations involving Google Cloud Storage!
-
-        To continue, please set the 'GOOGLE_APPLICATION_CREDENTIALS' variable
-        either:
-
-        1. In your terminal:
-        export GOOGLE_APPLICATION_CREDENTIALS='path-to-your-credentials-file.json'
-
-        2. In python:
-        import os
-        os.environ['GOOGLE_APPLICATION_CREDENTIALS'] = 'path-to-your-credentials-file.json'
-
-        Then import io once this is set
-        """
-    )
-
-## cannot get project automatically, so it is provided as argument
-storage_client = storage.Client(project = 'goes-olr')
-goes_16_bucket = storage_client.get_bucket("gcp-public-data-goes-16")
-
-
-def _test_subprocess_command(shell_command):
-    """
-    Test if a shell command can be successfully executed by suprocess
-
-    inputs:
-    -- shell_command (str): shell command that will be passed to subprocess
-
-    outputs:
-    -- boolean: True if command executed without errors. False if called process
-        returned an error
-    """
-    try:
-        test = subprocess.check_output(shell_command.split(" "))
-    except subprocess.CalledProcessError:
-        return False
-    else:
-        return True
-
-
-# Check if we can use ncdump to check if files are correct:
-ncdump_exists = _test_subprocess_command("ncdump")
-
-if not ncdump_exists:
-    warnings.warn(
-        """Warning: ncdump shell command not found, resorting to xarray
-                     for file checking"""
-    )
-
-
-def _check_file_size_against_blob(filename, blob, rtol=0.1, atol=1):
-    """
-    Check if a file size on disk is within toterance of a GCS blob
-
-    inputs:
-    -- filename (str): path to file
-    -- blob (Blob): GCS blob of file to compare
-    -- rtol (float; default=0.1): Relative tolerance in bytes (see numpy.isclose)
-    -- atol (float; default=1): Absolute tolerance in bytes (see numpy.isclose)
-
-    outputs:
-    -- boolean: True if file size is within tolerance of the blob. False if
-        otherwise
-    """
-    filesize = os.stat(filename).st_size
-    blobsize = blob.size
-    return np.isclose(filesize, blobsize, rtol, atol)
-
-
-def _check_ncdump_is_valid(filename):
-    """
-    Checks if a netcdf file can be examined using ncdump
-
-    inputs:
-    -- filename (str): path to file
-
-    outputs:
-    -- boolean: True if file can be examined using ncdump. False if otherwise
-    """
-    if ncdump_exists:
-        if _test_subprocess_command(f"ncdump -h {filename}"):
-            return True
-        else:
-            return False
-    else:
-        raise RuntimeError("ncdump cannot be called")
-
-
-def _check_xarray_is_valid(filename):
-    """
-    Checks if a netcdf file can be opened using xarray
-
-    inputs:
-    -- filename (str): path to file
-
-    outputs:
-    -- boolean: True if file can be opened using xarray. False if otherwise
-    """
-    try:
-        with xr.open_dataset(filename) as ds:
-            pass
-    except OSError:
-        return False
-    else:
-        return True
-
-
-def _check_netcdf_file_is_valid(filename):
-    """
-    Checks if a netcdf file can be opened. Uses ncdump if available, otherwise
-    attempts to open the file using xarrays
-
-    inputs:
-    -- filename (str): path to file
-
-    outputs:
-    -- boolean: True if file can be opened. False if otherwise
-    """
-    if ncdump_exists:
-        return _check_ncdump_is_valid(filename)
-    else:
-        return _check_xarray_is_valid(filename)
-
-
-def _check_free_space_for_blob(
-    blob, save_path, relative_min_storage=1.25, absolute_min_storage=2**30
-):
-    """
-    Checks if the avaialbe space on a storage volume is large enough to download
-    a blob. Min storage required is the largest out of the relative_min_storage
-    multiplied by the blob size, and the absolute_min_storage.
-
-    inputs:
-    -- blob (Blob): GCS blob of file to compare
-    -- save_path (str): direfctory to save file to
-    -- relative_min_storage (float; default=1.25): Relative tolerance in bytes (see numpy.isclose)
-    -- absolute_min_storage (float; default=2**30 (1GB)): Absolute tolerance in bytes (see numpy.isclose)
-
-    outputs:
-    -- boolean: True if free storage is above min threhsold, otherwise False.
-    """
-    min_storage = np.maximum(blob.size * relative_min_storage, absolute_min_storage)
-
-    if shutil.disk_usage(save_path).free > min_storage:
-        return True
-    else:
-        return False
-
-
-def _check_if_file_exists_and_is_valid(filename, blob, remove_corrupt=True):
-    """
-    Checks if a netcdf file exists and can be opened.
-
-    inputs:
-    -- filename (str): path to file
-    -- blob (Blob): GCS blob of file to compare
-    -- remove_corrupt (bool; default is True): If true, if the file exists but
-        is not valid then delete the file. Otherwise leave in place
-
-    outputs:
-    -- boolean: True if file can be opened. False if otherwise
-    """
-    if os.path.exists(filename):
-        if _check_file_size_against_blob(filename, blob):
-            return True
-        else:
-            if remove_corrupt:
-                os.remove(filename)
-
-            return False
-    else:
-        return False
-
-
-def _find_abi_blobs(date, satellite=16, product="Rad", view="C", mode=3, channel=1):
-    """
-    Find ABI file blobs for a single set of input parameters
-
-    inputs:
-    -- date (datetime): date to find files (hour)
-    -- satellite (int; default=16): GOES satellite to get data for. 16 or 17
-    -- product (str; default='Rad'): ABI data product to download
-    -- view (str; default='C'): View to get data for
-    -- mode (int; default=4): Scanning mode to get data for
-    -- channel (int; default=1): Channel to get data for. Only used for Rad and
-        CMIP products
-
-    outputs:
-    -- list of Blobs: list of blobs found using the prefix generated from the
-        supplied inputs
-    """
-    if satellite == 16:
-        goes_bucket = goes_16_bucket
-    elif satellite == 17:
-        goes_bucket = goes_17_bucket
-    else:
-        raise ValueError("Invalid input for satellite keyword")
-
-    doy = (date - datetime(date.year, 1, 1)).days + 1
-
-    level = "L1b" if product == "Rad" else "L2"
-    
-    blob_path = "ABI-%s-%s%.1s/%04d/%03d/%02d/" % (
-        level,
-        product,
-        view,
-        date.year,
-        doy,
-        date.hour,
-    )
-
-    if product == "Rad" or "CMIP":
-        blob_prefix = "OR_ABI-%s-%s%s-M%1dC%02d_G%2d_s" % (
-            level,
-            product,
-            view,
-            mode,
-            channel,
-            satellite,
-        )
-        
-    else:
-        blob_prefix = "OR_ABI-%s-%s%s-M%1d_G%2d_s" % (
-            level,
-            product,
-            view,
-            mode,
-            satellite,
-        )
-    print(f'path: {blob_path + blob_prefix}\n')
-    blobs = list(goes_bucket.list_blobs(prefix=blob_path + blob_prefix, delimiter="/"))
-
-    print(f'found blobs {blobs}\n\n')
-    print(f'returning first blob {[blobs[0]]}\n\n')
-    # return blobs
-    return [blobs[0]]
-
-
-def find_abi_blobs(
-    dates, satellite=16, product="Rad", view="C", mode=[3, 4, 6], channel=1
-):
-    """
-    Find ABI file blobs for a set of one or multiple input parameters
-
-    inputs:
-    -- dates (datetime): date or list of dates to find files (hour)
-    -- satellite (int; default=16): GOES satellite to get data for. 16 or 17
-    -- product (str; default='Rad'): ABI data product to download
-    -- view (str; default='C'): View to get data for
-    -- mode (int; default=[3,4,6]): Scanning mode to get data for
-    -- channel (int; default=1): Channel to get data for. Only used for Rad and
-        CMIP products
-
-    outputs:
-    -- list of Blobs: list of blobs found using the prefix generated from the
-        supplied inputs
-    """
-
-    input_params = [
-        m.ravel().tolist()
-        for m in np.meshgrid(dates, satellite, product, view, mode, channel)
-    ]
-    n_params = len(input_params[0])
-
-    blobs = sum(
-        [
-            _find_abi_blobs(
-                input_params[0][i],
-                satellite=input_params[1][i],
-                product=input_params[2][i],
-                view=input_params[3][i],
-                mode=input_params[4][i],
-                channel=input_params[5][i],
-            )
-            for i in range(n_params)
-        ],
-        [],
-    )
-
-    return blobs
-
-
-def _get_download_destination(blob, save_dir, replicate_path=True):
-    """
-    Generate path to download destination for a blob
-
-    inputs:
-    -- blob (Blob): GCS blob to find download location for
-    -- save_dir (str): root directory to save files in
-    -- replicate_path (bool; default=True): If True, replicate the directory
-        structure within the GCS store. Otherwise, provide a path to the file
-        directly within the root directory
-
-    outputs:
-    -- str: save_file name and location for download
-    """
-    blob_path, blob_name = os.path.split(blob.name)
-
-    if replicate_path:
-        save_path = os.path.join(save_dir, blob_path)
-    else:
-        save_path = save_dir
-    if not os.path.isdir(save_path):
-        os.makedirs(save_path, exist_ok=True)
-    save_file = os.path.join(save_path, blob_name)
-    return save_file
-
-
-def download_blob(
-    blob,
-    save_dir,
-    replicate_path=True,
-    check_download=False,
-    n_attempts=1,
-    clobber=False,
-    min_storage=2**30,
-    verbose=False,
-    remove_corrupt=True,
-):
-    """
-    Download a single blob from GCS
-
-    inputs:
-    -- blob (Blob): GCS blob to find download location for
-    -- save_dir (str): root directory to save files in
-    -- replicate_path (bool; default=True): If True, replicate the directory
-        structure within the GCS store. Otherwise, provide a path to the file
-        directly within the root directory
-    -- check_download (bool; default=True)
-    """
-    save_path = _get_download_destination(blob, save_dir, replicate_path=replicate_path)
-    if clobber or not os.path.exists(save_path):
-        if clobber and os.path.exists(save_path):
-            os.remove(save_path)
-        if verbose:
-            print(f"Downloading {save_path}", flush=True)
-        if _check_free_space_for_blob(
-            blob,
-            os.path.split(save_path)[0],
-            relative_min_storage=1.25,
-            absolute_min_storage=min_storage,
-        ):
-            try:
-                blob.download_to_filename(save_path)
-            except OSError:
-                if n_attempts > 1:
-                    download_blob(
-                        blob,
-                        save_dir,
-                        replicate_path=replicate_path,
-                        check_download=check_download,
-                        n_attempts=n_attempts - 1,
-                        clobber=clobber,
-                    )
-                else:
-                    raise RuntimeError(f"{save_path}: download failed")
-            if check_download and not _check_if_file_exists_and_is_valid(
-                save_path, blob, remove_corrupt=True
-            ):
-                if n_attempts > 1:
-                    download_blob(
-                        blob,
-                        save_dir,
-                        replicate_path=replicate_path,
-                        check_download=check_download,
-                        n_attempts=n_attempts - 1,
-                        clobber=clobber,
-                    )
-                else:
-                    if remove_corrupt:
-                        if os.path.exists(save_path):
-                            os.remove(save_path)
-                    raise RuntimeError(f"{save_path}: downloaded file not valid")
-        else:
-            raise OSError("Not enough storage space available for download")
-    if os.path.exists(save_path):
-        return save_path
-    else:
-        raise RuntimeError(f"{save_path}: downloaded file not found")
-
-
-
-def get_goes_date(filename: str) -> datetime:
-    """
-    Finds the centre point time from an ABI filename
-    """
-    base_string = os.path.split(filename)[-1]
-
-    start_string = base_string.split("_s")[-1]
-    start_date = parse_date(start_string[:4] + "0101" + start_string[7:13]) + timedelta(
-        days=int(start_string[4:7]) - 1
-    )
-    end_string = base_string.split("_e")[-1]
-    end_date = parse_date(end_string[:4] + "0101" + end_string[7:13]) + timedelta(
-        days=int(end_string[4:7]) - 1
-    )
-
-    return start_date + (end_date - start_date) / 2
-
-
-def find_abi_files(
-    date,
-    satellite=16,
-    product="Rad",
-    view="C",
-    mode=[3, 4, 6],
-    channel=1,
-    save_dir="./",
-    replicate_path=True,
-    check_download=False,
-    n_attempts=1,
-    download_missing=False,
-    clobber=False,
-    min_storage=2**30,
-    remove_corrupt=True,
-    verbose=False,
-):
-    """
-    Finds ABI files on the local system. Optionally downloads missing files from
-    GCS
-
-    inputs:
-    # TODO:
-
-    outputs:
-        list: list of filenames on local system
-    """
-    blobs = find_abi_blobs(
-        date,
-        satellite=satellite,
-        product=product,
-        view=view,
-        mode=mode,
-        channel=channel,
-    )
-    files = []
-    for blob in blobs:
-        if download_missing:
-            try:
-                save_file = download_blob(
-                    blob,
-                    save_dir,
-                    replicate_path=replicate_path,
-                    check_download=check_download,
-                    n_attempts=n_attempts,
-                    clobber=clobber,
-                    min_storage=min_storage,
-                    verbose=verbose,
-                )
-            except OSError as e:
-                warnings.warn(str(e.args[0]))
-                download_missing = False
-            except RuntimeError as e:
-                warnings.warn(str(e.args[0]))
-            else:
-                if os.path.exists(save_file):
-                    files += [save_file]
-        else:
-            local_file = _get_download_destination(
-                blob, save_dir, replicate_path=replicate_path
-            )
-            if _check_if_file_exists_and_is_valid(
-                local_file, blob, remove_corrupt=remove_corrupt
-            ):
-                files += [local_file]
-    return files
-
-
-def _find_glm_blobs(date, satellite=16):
-    if satellite == 16:
-        goes_bucket = goes_16_bucket
-    elif satellite == 17:
-        goes_bucket = goes_17_bucket
-    else:
-        raise ValueError("Invalid input for satellite keyword")
-
-    doy = (date - datetime(date.year, 1, 1)).days + 1
-
-    blob_path = "GLM-L2-LCFA/%04d/%03d/%02d/" % (date.year, doy, date.hour)
-    blob_prefix = "OR_GLM-L2-LCFA_G%2d_s" % satellite
-
-    blobs = list(goes_bucket.list_blobs(prefix=blob_path + blob_prefix, delimiter="/"))
-
-    return blobs
-
-
-def find_glm_blobs(dates, satellite=16):
-    """
-    Find GLM file blobs for a set of one or multiple input parameters
-
-    inputs:
-    -- dates (datetime): date or list of dates to find files (hour)
-    -- satellite (int; default=16): GOES satellite to get data for. 16 or 17
-
-    outputs:
-    -- list of Blobs: list of blobs found using the prefix generated from the
-        supplied inputs
-    """
-
-    input_params = [m.ravel().tolist() for m in np.meshgrid(dates, satellite)]
-    n_params = len(input_params[0])
-
-    blobs = sum(
-        [
-            _find_glm_blobs(input_params[0][i], satellite=input_params[1][i])
-            for i in range(n_params)
-        ],
-        [],
-    )
-
-    return blobs
-
-
-def find_glm_files(
-    date,
-    satellite=16,
-    save_dir="./",
-    replicate_path=True,
-    check_download=False,
-    n_attempts=1,
-    download_missing=False,
-    clobber=False,
-    min_storage=2**30,
-    remove_corrupt=True,
-    verbose=False,
-):
-    """
-    Finds GLM files on the local system. Optionally downloads missing files from
-    GCS
-
-    inputs:
-    # TODO:
-
-    outputs:
-        list: list of filenames on local system
-    """
-    blobs = find_glm_blobs(date, satellite=satellite)
-    files = []
-    for blob in blobs:
-        if download_missing:
-            try:
-                save_file = download_blob(
-                    blob,
-                    save_dir,
-                    replicate_path=replicate_path,
-                    check_download=check_download,
-                    n_attempts=n_attempts,
-                    clobber=clobber,
-                    min_storage=min_storage,
-                    verbose=verbose,
-                )
-            except OSError as e:
-                warnings.warn(e.args[0])
-                download_missing = False
-            except RuntimeError:
-                warnings.warn(e.args[0])
-            else:
-                if os.path.exists(save_file):
-                    files += [save_file]
-        else:
-            save_file = _get_download_destination(
-                blob, save_dir, replicate_path=replicate_path
-            )
-            if _check_if_file_exists_and_is_valid(save_file, blob):
-                files += [save_file]
-    return files
-
-
-# def find_nexrad_blobs(date, site):
-#     storage_client = storage.Client()
-#     bucket = storage_client.get_bucket("gcp-public-data-nexrad-l2")
-
-#     blob_path = "%04d/%02d/%02d/%s/" % (date.year, date.month, date.day, site)
-#     blob_prefix = "NWS_NEXRAD_NXL2DPBL_%s_%04d%02d%02d%02d" % (
-#         site,
-#         date.year,
-#         date.month,
-#         date.day,
-#         date.hour,
-#     )
-
-#     blobs = list(bucket.list_blobs(prefix=blob_path + blob_prefix, delimiter="/"))
-
-#     return blobs
-
-
-def download_blobs(
-    blob_list, save_dir="./", replicate_path=True, n_attempts=0, clobber=False
-):
-    for blob in blob_list:
-        blob_path, blob_name = os.path.split(blob.name)
-
-        if replicate_path:
-            save_path = os.path.join(save_dir, blob_path)
-        else:
-            save_path = save_dir
-        if not os.path.isdir(save_path):
-            os.makedirs(save_path, exist_ok=True)
-
-        save_file = os.path.join(save_path, blob_name)
-        if clobber or not os.path.exists(save_file):
-            blob.download_to_filename(save_file)
-
-
-# def find_nexrad_files(
-#     date, site, save_dir="./", replicate_path=True, download_missing=False
-# ):
-#     blobs = find_nexrad_blobs(date, site)
-#     files = []
-#     for blob in blobs:
-#         blob_path, blob_name = os.path.split(blob.name)
-
-#         if replicate_path:
-#             save_path = os.path.join(save_dir, blob_path)
-#         else:
-#             save_path = save_dir
-#         if not os.path.isdir(save_path):
-#             os.makedirs(save_path, exist_ok=True)
-
-#         save_file = os.path.join(save_path, blob_name)
-#         if os.path.exists(save_file):
-#             files += [save_file]
-#         elif download_missing:
-#             download_blobs([blob], save_dir=save_dir, replicate_path=replicate_path)
-#             if os.path.exists(save_file):
-#                 files += [save_file]
-
-#     return files
\ No newline at end of file
diff --git a/notebooks/goes-scripts/xy_to_latlon.py b/notebooks/goes-scripts/xy_to_latlon.py
deleted file mode 100644
index 806d155..0000000
--- a/notebooks/goes-scripts/xy_to_latlon.py
+++ /dev/null
@@ -1,73 +0,0 @@
-import numpy as np
-
-'''
-script adapted from https://lsterzinger.medium.com/add-lat-lon-coordinates-to-goes-16-goes-17-l2-data-and-plot-with-cartopy-27f07879157f
-'''
-
-
-def calc_latlon(ds):
-    '''
-    Takes GOES dataset (one image) and computes latitude and 
-    longitude for each pixel using horizontal scan angles x
-    and vertical scan angles y.
-
-    Input:
-        ds  xarray.Dataset
-    Output:
-        ds  xarray.Dataset with lat and lon values added for each datapoint
-            and used as indeces.
-    '''
-    # The math for this function was taken from 
-    # https://makersportal.com/blog/2018/11/25/goes-r-satellite-latitude-and-longitude-grid-projection-algorithm 
-
-    x = ds.x
-    y = ds.y
-    goes_imager_projection = ds.goes_imager_projection
-    
-    x,y = np.meshgrid(x,y)
-    
-    r_eq = goes_imager_projection.attrs["semi_major_axis"] # earth radius at equator
-    r_pol = goes_imager_projection.attrs["semi_minor_axis"] # earth radius at pole
-    l_0 = goes_imager_projection.attrs["longitude_of_projection_origin"] * (np.pi/180) # lambda0
-    h_sat = goes_imager_projection.attrs["perspective_point_height"] # distance satellite to nearest equator surface point
-    H = r_eq + h_sat # distance satellite to earth centre
-    
-    a = np.sin(x)**2 + (np.cos(x)**2 * (np.cos(y)**2 + (r_eq**2 / r_pol**2) * np.sin(y)**2))
-    b = -2 * H * np.cos(x) * np.cos(y)
-    c = H**2 - r_eq**2
-    
-    r_s = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)
-    
-    s_x = r_s * np.cos(x) * np.cos(y)
-    s_y = -r_s * np.sin(x)
-    s_z = r_s * np.cos(x) * np.sin(y)
-
-    # latitude and longitude
-    lat = np.arctan((r_eq**2 / r_pol**2) * (s_z / np.sqrt((H-s_x)**2 +s_y**2))) * (180/np.pi)
-    lon = (l_0 - np.arctan(s_y / (H-s_x))) * (180/np.pi)
-    
-    ds = ds.assign_coords({
-        "lat":(["y","x"],lat),
-        "lon":(["y","x"],lon)
-    })
-    ds.lat.attrs["units"] = "degrees_north"
-    ds.lon.attrs["units"] = "degrees_east"
-
-    return ds
-
-def get_xy_from_latlon(ds, lats, lons):
-    lat1, lat2 = lats
-    lon1, lon2 = lons
-
-    lat = ds.lat.data
-    lon = ds.lon.data
-    
-    x = ds.x.data
-    y = ds.y.data
-    
-    x,y = np.meshgrid(x,y)
-    
-    x = x[(lat >= lat1) & (lat <= lat2) & (lon >= lon1) & (lon <= lon2)]
-    y = y[(lat >= lat1) & (lat <= lat2) & (lon >= lon1) & (lon <= lon2)] 
-    
-    return ((min(x), max(x)), (min(y), max(y)))
diff --git a/rs_tools/_src/geoprocessing/grid.py b/rs_tools/_src/geoprocessing/grid.py
index 4ca30ba..70b90d4 100644
--- a/rs_tools/_src/geoprocessing/grid.py
+++ b/rs_tools/_src/geoprocessing/grid.py
@@ -1,7 +1,7 @@
 import numpy as np
 from typing import Tuple
 
-
+# TODO: To be moved to Earth System Datacube Tools
 def create_latlon_grid(region: Tuple[float, float, float, float],
                        resolution: float):
     """
diff --git a/rs_tools/_src/geoprocessing/modis/interp.py b/rs_tools/_src/geoprocessing/modis/interp.py
index 1fa8ce1..8e1bf9c 100644
--- a/rs_tools/_src/geoprocessing/modis/interp.py
+++ b/rs_tools/_src/geoprocessing/modis/interp.py
@@ -2,6 +2,7 @@
 import numpy as np
 from scipy.interpolate import griddata
 
+# TODO: Satpy can interpolate grid, so function not necessarily needed
 def interp_coords_modis(ds: xr.Dataset, desired_size: tuple[int], method: str="cubic") -> np.array:
 
     lat = ds.Latitude
diff --git a/rs_tools/_src/geoprocessing/modis/rescale.py b/rs_tools/_src/geoprocessing/modis/rescale.py
index b9b95b8..93b0952 100644
--- a/rs_tools/_src/geoprocessing/modis/rescale.py
+++ b/rs_tools/_src/geoprocessing/modis/rescale.py
@@ -2,6 +2,7 @@
 import xarray as xr
 import numpy as np
 
+# TODO: Satpy can rescale data, so function not necessarily needed
 def convert_integers2radiances(
         da: xr.DataArray, 
     ) -> xr.DataArray:
diff --git a/rs_tools/_src/geoprocessing/reproject.py b/rs_tools/_src/geoprocessing/reproject.py
index abd4e1e..82c3ce9 100644
--- a/rs_tools/_src/geoprocessing/reproject.py
+++ b/rs_tools/_src/geoprocessing/reproject.py
@@ -3,6 +3,7 @@
 from pyproj import CRS, Transformer
 import xarray as xr
 
+# TODO: To be moved to Earth System Datacube Tools
 def convert_lat_lon_to_x_y(crs: str, lon: list[float], lat: list[float]) -> Tuple[float, float]:
     transformer = Transformer.from_crs("epsg:4326", crs, always_xy=True)
     x, y = transformer.transform(lon, lat)
diff --git a/rs_tools/_src/modules/__init__.py b/rs_tools/_src/modules/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/rs_tools/_src/modules/modis/__init__.py b/rs_tools/_src/modules/modis/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/rs_tools/_src/modules/modis/download.py b/rs_tools/_src/modules/modis/download.py
deleted file mode 100644
index e69de29..0000000
diff --git a/rs_tools/_src/preprocessing/normalize.py b/rs_tools/_src/preprocessing/normalize.py
index 9a46f57..0e618c9 100644
--- a/rs_tools/_src/preprocessing/normalize.py
+++ b/rs_tools/_src/preprocessing/normalize.py
@@ -7,11 +7,11 @@
 def spatial_mean(ds: xr.Dataset, spatial_variables: List[str]) -> xr.Dataset:
     return ds.mean(spatial_variables)
 
-
+# TODO: Check this function
 def normalize(
-        files: List[str], 
+        files: List[str],
+        temporal_variables: List[str], 
         spatial_variables: List[str]=["x","y"], 
-        temporal_variables: List[str],
 ) -> xr.Dataset:
     
     preprocess = partial(spatial_mean, spatial_variables=spatial_variables)
diff --git a/rs_tools/_src/utils/crs_ops.py b/rs_tools/_src/utils/crs_ops.py
deleted file mode 100644
index bc7f2ab..0000000
--- a/rs_tools/_src/utils/crs_ops.py
+++ /dev/null
@@ -1,120 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Oct 20 23:49:15 2020
-
-@author: ghiggi
-"""
-import numpy as np
-import cartopy.crs as ccrs
-
-def get_bbox_polygon(bbox): 
-    from shapely.geometry.polygon import Polygon
-    p = Polygon(((bbox[0],bbox[1]), # bottom left
-                (bbox[0],bbox[3]),  # top left
-                (bbox[2],bbox[3]),  # top rght
-                (bbox[2],bbox[1]),  # bottom right
-                (bbox[0],bbox[1]))) # bottom left
-    return p
-
-def get_bbox_coords(bbox, output_type='xy'):
-    """Return corner coordinates of bounding box.
-    # bbox : (xmin, ymin, xmax, ymax) 
-    
-    # output_type = 'shapely', 'array', 'xy', 'tuple'
-    # shapely: shapely polygon 
-    # array : numpy array with xy columns
-    # xy : x, y numpy arrays 
-    # tuple: list of (x,y)  
-    
-    # Return (bottom_left, top_left, top_right bottom_right)
-    """
-    if (output_type not in ['xy','tuple','array']):
-        raise ValueError("Valid output_type: 'xy','tuple','array'")
-    x = np.array([bbox[0], bbox[0], bbox[2], bbox[2]])
-    y = np.array([bbox[1], bbox[3], bbox[3], bbox[1]])
-    if (output_type=='xy'):
-        return x, y
-    if (output_type=='tuple'):
-        return list(zip(x,y))
-    if (output_type=='array'):
-        return np.vstack((x,y)).transpose()
-    
-
-def get_bbox(x,y):
-    xmin = np.min(x)
-    xmax = np.max(x)
-    ymin = np.min(y)
-    ymax = np.max(y)
-    return [xmin, ymin, xmax, ymax]
-
-def get_projected_bbox(bbox, crs_ref, crs_proj, lonlat=False): 
-    # Check if xmin > xmax (crossing the antimeridian)
-    # - This only occurs for geographic coordinates / projections
-    # TODO : Infer if lonlat from CRS 
-    flag_crossing = False 
-    if (bbox[0] > bbox[2]):
-        flag_crossing = True
-    #------------------------------------------------------------.   
-    x, y = get_bbox_coords(bbox, output_type="xy")  
-    xy_prj = crs_proj.transform_points(src_crs=crs_ref, x=x, y=y)[:,0:2]
-    bbox_proj = get_bbox(xy_prj[:,0], xy_prj[:,1]) # this simply look at min and max
-    #------------------------------------------------------------.  
-    # Deal with crossing antimeridian
-    if (flag_crossing is True) and (lonlat is True): 
-        xmax = bbox_proj[0]
-        xmin = bbox_proj[2]
-        bbox_proj[0] = xmin 
-        bbox_proj[2] = xmax 
-    #------------------------------------------------------------.      
-    return bbox_proj
-
-
-def adapt_bbox_for_resolution(bbox, resolution, lonlat=False):
-    """ Adapt bbox to allow specified resolution. """
-    # Check resolution input 
-    if (not isinstance(resolution, (float, int, tuple))):
-        raise TypeError("Provide resolution as single integer/float or as tuple.")
-    if (isinstance(resolution, (float, int))):
-        resolution = (resolution,resolution)
-    else:
-        if len(resolution) != 2:
-            raise ValueError("Please provide resolution in (x,y) tuple format.")
-    # Retrieve distance in x and y direction 
-    dx = bbox[2] - bbox[0]
-    dy = bbox[3] - bbox[1]
-    # Check dx and dy are larger than resolution 
-    if (dx < resolution[0]):
-        raise ValueError("The specified resolution is larger than the bbox in x direction")
-    if (dy < resolution[1]):
-        raise ValueError("The specified resolution is larger than the bbox in y direction")    
-    #-------------------------------------------------------------------------.
-    # Retrieve corrections
-    x_corr = (np.ceil(dx/resolution[0])*resolution[0] - dx)/2
-    y_corr = (np.ceil(dy/resolution[1])*resolution[1] - dy)/2
-    xmin = bbox[0]-x_corr
-    ymin = bbox[1]-y_corr
-    xmax = bbox[2]+x_corr
-    ymax = bbox[3]+y_corr
-    #-------------------------------------------------------------------------.
-    # Check y is outside [-90, 90] in lat ...
-    if (lonlat is True):
-        if (ymin < -90): 
-            thr = -90-ymin
-            ymin = -90 
-            ymax = ymax + thr
-        if (ymax > 90): 
-            thr = ymax -90 
-            ymax = 90 
-            ymin = ymin - thr
-            if (ymin > -90):
-                print("Impossible to adapt the bbox by extending it.")
-                ymin = -90
-    #-------------------------------------------------------------------------. 
-    # Ensure that xmin and xmax stay inside [-180, 180] boundary
-    if (xmin > 180): xmin = xmin - 360
-    if (xmax > 180): xmax = xmax - 360
-    if (xmin < 180): xmin = xmin + 360
-    if (xmax < 180): xmax = xmax + 360
-    #-------------------------------------------------------------------------. 
-    return [xmin,ymin, xmax, ymax]
diff --git a/rs_tools/_src/utils/math.py b/rs_tools/_src/utils/math.py
index c774567..e16a002 100644
--- a/rs_tools/_src/utils/math.py
+++ b/rs_tools/_src/utils/math.py
@@ -1,5 +1,6 @@
 from math import floor, ceil
 
+# TODO: Eman to refactor code
 def bounds_and_points_to_step(xmin: float, xmax: float, Nx: float) -> float:
     """Calculates the dx from the minmax
     Eq: