diff --git a/pylst/__init__.py b/pylst/__init__.py
index 4dc3c21..2285a57 100644
--- a/pylst/__init__.py
+++ b/pylst/__init__.py
@@ -1,5 +1,37 @@
-"""Top-level package for pylst."""
+# __init__.py
 
-__author__ = """Azad Rasul"""
-__email__ = 'azad.rasul@soran.edu.iq'
-__version__ = '0.0.1'
+# Importing specific functionalities from different modules within the package
+
+# Importing from the 'pylst' module
+from .pylst import emissivity
+
+# Importing visualization functions from the 'visualization' module
+from .visualization import (
+    plot_data,        # Function to plot data
+    show,             # Function to show visualization
+    rastshow,         # Function to show raster visualization
+    create_interactive_map,  # Function to create interactive map
+    histogram_equalization   # Function for histogram equalization
+)
+
+# Importing spatial analysis functions from the 'spatial_analysis' module
+from .spatial_analysis.zonstat import calculate_zs  # Zonal statistics
+from .spatial_analysis.changedet import analyze_images  # Change Detection
+
+# Importing image handling classes and functions from the 'open' module
+from .open import (
+    Opener,        # Class for opening image data
+    remopener,     # Class for opening remote image data
+    download_images  # Function to download image collection images
+)
+
+# Importing normalization function from the 'normalization' module
+from .normalization import nrs  # Normalization Ratio Scale function
+
+# Importing machine learning functions from the 'ml' module
+from .ml import (
+    train_regression_model,      # Function to train regression model
+    train_classification_model,  # Function to train classification model
+    train_test_split,            # Function to split data for training and testing
+    accuracy_score               # Function to compute accuracy score
+)
diff --git a/pylst/emissivity/__init__.py b/pylst/emissivity/__init__.py
new file mode 100644
index 0000000..1cd09c9
--- /dev/null
+++ b/pylst/emissivity/__init__.py
@@ -0,0 +1,2 @@
+# Importing the default algorithms for emissivity calculations
+from .emissivity import default_algorithms
diff --git a/pylst/emissivity/algorithms.py b/pylst/emissivity/algorithms.py
new file mode 100644
index 0000000..328b2fa
--- /dev/null
+++ b/pylst/emissivity/algorithms.py
@@ -0,0 +1,242 @@
+import numpy as np
+
+from pylandtemp.utils import rescale_band, cavity_effect, fractional_vegetation_cover
+
+# Credit to the original author: Oladimeji Mudele (pylandtemp)
+# Original source: https://github.com/pylandtemp/pylandtemp
+
+class Emissivity:
+    def __init__(self):
+        """Parent class for all emissivity methods. Contains general methods and attributes"""
+        self.ndvi_min = -1
+        self.ndvi_max = 1
+        self.baresoil_ndvi_max = 0.2
+        self.vegatation_ndvi_min = 0.5
+
+    def __call__(self, **kwargs) -> np.ndarray:
+        """Computes the emissivity
+
+        kwargs:
+        **ndvi (np.ndarray): NDVI image
+        **red_band (np.ndarray): Band 4 or Red band image.
+        **mask (np.ndarray[bool]): Mask image. Output will have NaN value where mask is True.
+
+
+        Returns:
+            Tuple(np.ndarray, np.ndarray): Emissivity for bands 10 and 11 respectively
+        """
+        if "ndvi" not in kwargs:
+            raise ValueError("NDVI image is not provided")
+
+        if "red_band" not in kwargs:
+            raise ValueError("Band 4 (red band) image is not provided")
+
+        self.ndvi = kwargs["ndvi"]
+        self.red_band = kwargs["red_band"]
+
+        if (
+            self.ndvi is not None
+            and self.red_band is not None
+            and self.ndvi.shape != self.red_band.shape
+        ):
+            raise ValueError(
+                "Input images (NDVI and Red band) must be of equal dimension"
+            )
+
+        emm_10, emm_11 = self._compute_emissivity()
+        mask = emm_10 == 0
+        emm_10[mask] = np.nan
+        if emm_11 is not None:
+            emm_11[mask] = np.nan
+        return emm_10, emm_11
+
+    def _compute_emissivity(self):
+        """Abstract method to be implemented by subclasses"""
+        raise NotImplementedError("No concrete implementation of emissivity method yet")
+
+    def _get_land_surface_mask(self):
+        """Return a mask for differnet land surface types"""
+        mask_baresoil = (self.ndvi >= self.ndvi_min) & (
+            self.ndvi < self.baresoil_ndvi_max
+        )
+        mask_vegetation = (self.ndvi > self.vegatation_ndvi_min) & (
+            self.ndvi <= self.ndvi_max
+        )
+        mask_mixed = (self.ndvi >= self.baresoil_ndvi_max) & (
+            self.ndvi <= self.vegatation_ndvi_min
+        )
+        return {
+            "baresoil": mask_baresoil,
+            "vegetation": mask_vegetation,
+            "mixed": mask_mixed,
+        }
+
+    def _get_landcover_mask_indices(self):
+        """Returns indices corresponding to the different landcover classes of of interest namely:
+        vegetation, baresoil and mixed"
+        """
+        masks = self._get_land_surface_mask()
+        baresoil = np.where(masks["baresoil"])
+        vegetation = np.where(masks["vegetation"])
+        mixed = np.where(masks["mixed"])
+        return {"baresoil": baresoil, "vegetation": vegetation, "mixed": mixed}
+
+    def _compute_fvc(self):
+        # Returns the fractional vegegation cover from the NDVI image.
+        return fractional_vegetation_cover(self.ndvi)
+
+
+class ComputeMonoWindowEmissivity(Emissivity):
+
+    emissivity_soil_10 = 0.97
+    emissivity_veg_10 = 0.99
+    emissivity_soil_11 = None
+    emissivity_veg_11 = None
+
+    def _compute_emissivity(self) -> np.ndarray:
+        emm = np.empty_like(self.ndvi)
+        landcover_mask_indices = self._get_landcover_mask_indices()
+
+        # Baresoil value assignment
+        emm[landcover_mask_indices["baresoil"]] = self.emissivity_soil_10
+        # Vegetation value assignment
+        emm[landcover_mask_indices["vegetation"]] = self.emissivity_veg_10
+        # Mixed value assignment
+        emm[landcover_mask_indices["mixed"]] = (
+            0.004
+            * (((self.ndvi[landcover_mask_indices["mixed"]] - 0.2) / (0.5 - 0.2)) ** 2)
+        ) + 0.986
+        return emm, emm
+
+
+class ComputeEmissivityNBEM(Emissivity):
+    """
+    Method references:
+
+    1. Li, Tianyu, and Qingmin Meng. "A mixture emissivity analysis method for
+        urban land surface temperature retrieval from Landsat 8 data." Landscape
+        and Urban Planning 179 (2018): 63-71.
+
+    2. Yu, Xiaolei, Xulin Guo, and Zhaocong Wu. "Land surface temperature retrieval
+        from Landsat 8 TIRS—Comparison between radiative transfer equation-based method,
+        split window algorithm and single channel method." Remote sensing 6.10 (2014): 9829-9852.
+
+    """
+
+    emissivity_soil_10 = 0.9668
+    emissivity_veg_10 = 0.9863
+    emissivity_soil_11 = 0.9747
+    emissivity_veg_11 = 0.9896
+
+    def _compute_emissivity(self) -> np.ndarray:
+
+        if self.red_band is None:
+            raise ValueError(
+                "Red band cannot be {} for this emissivity computation method".format(
+                    self.red_band
+                )
+            )
+
+        self.red_band = rescale_band(self.red_band)
+        landcover_mask_indices = self._get_landcover_mask_indices()
+        fractional_veg_cover = self._compute_fvc()
+
+        def calc_emissivity_for_band(
+            image,
+            emissivity_veg,
+            emissivity_soil,
+            cavity_effect,
+            red_band_coeff_a=None,
+            red_band_coeff_b=None,
+        ):
+            image[landcover_mask_indices["baresoil"]] = red_band_coeff_a - (
+                red_band_coeff_b * self.red_band[landcover_mask_indices["baresoil"]]
+            )
+
+            image[landcover_mask_indices["mixed"]] = (
+                (emissivity_veg * fractional_veg_cover[landcover_mask_indices["mixed"]])
+                + (
+                    emissivity_soil
+                    * (1 - fractional_veg_cover[landcover_mask_indices["mixed"]])
+                )
+                + cavity_effect[landcover_mask_indices["mixed"]]
+            )
+
+            image[landcover_mask_indices["vegetation"]] = (
+                emissivity_veg + cavity_effect[landcover_mask_indices["vegetation"]]
+            )
+            return image
+
+        emissivity_band_10 = np.empty_like(self.ndvi)
+        emissivity_band_11 = np.empty_like(self.ndvi)
+        frac_vegetation_cover = self._compute_fvc()
+
+        cavity_effect_10 = cavity_effect(
+            self.emissivity_veg_10, self.emissivity_soil_10, fractional_veg_cover
+        )
+        cavity_effect_11 = cavity_effect(
+            self.emissivity_veg_11, self.emissivity_soil_11, fractional_veg_cover
+        )
+
+        emissivity_band_10 = calc_emissivity_for_band(
+            emissivity_band_10,
+            self.emissivity_veg_10,
+            self.emissivity_soil_10,
+            cavity_effect_10,
+            red_band_coeff_a=0.973,
+            red_band_coeff_b=0.047,
+        )
+        emissivity_band_11 = calc_emissivity_for_band(
+            emissivity_band_11,
+            self.emissivity_veg_11,
+            self.emissivity_soil_11,
+            cavity_effect_11,
+            red_band_coeff_a=0.984,
+            red_band_coeff_b=0.026,
+        )
+        return emissivity_band_10, emissivity_band_11
+
+
+class ComputeEmissivityGopinadh(Emissivity):
+    """
+    Method reference:
+
+    Rongali, Gopinadh, et al. "Split-window algorithm for retrieval of land surface temperature
+    using Landsat 8 thermal infrared data." Journal of Geovisualization and Spatial Analysis 2.2
+    (2018): 1-19.
+    """
+
+    emissivity_soil_10 = 0.971
+    emissivity_veg_10 = 0.987
+
+    emissivity_soil_11 = 0.977
+    emissivity_veg_11 = 0.989
+
+    def _compute_emissivity(self) -> np.ndarray:
+
+        fractional_veg_cover = self._compute_fvc()
+
+        def calc_emissivity_for_band(
+            image, emissivity_veg, emissivity_soil, fractional_veg_cover
+        ):
+            emm = (emissivity_soil * (1 - fractional_veg_cover)) + (
+                emissivity_veg * fractional_veg_cover
+            )
+            return emm
+
+        emissivity_band_10 = np.empty_like(self.ndvi)
+        emissivity_band_10 = calc_emissivity_for_band(
+            emissivity_band_10,
+            self.emissivity_veg_10,
+            self.emissivity_soil_10,
+            fractional_veg_cover,
+        )
+
+        emissivity_band_11 = np.empty_like(self.ndvi)
+        emissivity_band_11 = calc_emissivity_for_band(
+            emissivity_band_11,
+            self.emissivity_veg_11,
+            self.emissivity_soil_11,
+            fractional_veg_cover,
+        )
+        return emissivity_band_10, emissivity_band_11
\ No newline at end of file
diff --git a/pylst/emissivity/emissivity.py b/pylst/emissivity/emissivity.py
new file mode 100644
index 0000000..9307bd2
--- /dev/null
+++ b/pylst/emissivity/emissivity.py
@@ -0,0 +1,19 @@
+# Import emissivity computation algorithms from the algorithms module
+from .algorithms import (
+    ComputeMonoWindowEmissivity,
+    ComputeEmissivityNBEM,
+    ComputeEmissivityGopinadh,
+)
+
+# Define a dictionary of default emissivity computation algorithms
+# Each key corresponds to a method name, and the value is the corresponding algorithm class
+default_algorithms = {
+    "avdan": ComputeMonoWindowEmissivity,   # Avdan Ugur et al, 2016 method
+    "xiaolei": ComputeEmissivityNBEM,       # Xiaolei Yu et al, 2014 method
+    "gopinadh": ComputeEmissivityGopinadh,  # Gopinadh Rongali et al, 2018 method
+}
+
+# This dictionary allows users to easily choose a specific algorithm by providing the method name.
+# For example, to use the Avdan Ugur et al, 2016 method, users can select "avdan" as the method.
+# The corresponding algorithm class (ComputeMonoWindowEmissivity) will be used for computation.
+# Users can extend this dictionary by adding more methods and their corresponding algorithm classes.
diff --git a/pylst/exceptions.py b/pylst/exceptions.py
new file mode 100644
index 0000000..63910dc
--- /dev/null
+++ b/pylst/exceptions.py
@@ -0,0 +1,49 @@
+import numpy as np
+
+
+__all__ = [
+    "InvalidMaskError",
+    "InputShapesNotEqual",
+    "InvalidMethodRequested",
+]
+
+
+class InvalidMaskError(Exception):
+    """Invalid mask error"""
+
+    pass
+
+
+class KeywordArgumentError(Exception):
+    """Required keyword argument missing"""
+
+    pass
+
+
+class InputShapesNotEqual(Exception):
+    """Input images don't have the same shape"""
+
+    pass
+
+
+class InvalidMethodRequested(Exception):
+    """Invalid method/algorithm requested"""
+
+    pass
+
+
+def assert_required_keywords_provided(keywords, **kwargs):
+    """
+    This method checks if all the required keyword arguments to complete a computation
+        are provided in **kwargs
+    Args:
+        keywords ([list[str]], optional): Required keywords.
+    Raises:
+        KeywordArgumentError: custom exception
+    """
+    for keyword in keywords:
+        if keyword not in kwargs or kwargs[keyword] is None:
+            message = (
+                f"Keyword argument {keyword} must be provided for this computation "
+            )
+            raise KeywordArgumentError(message)
diff --git a/pylst/ml/__init__.py b/pylst/ml/__init__.py
new file mode 100644
index 0000000..6145fee
--- /dev/null
+++ b/pylst/ml/__init__.py
@@ -0,0 +1,7 @@
+# pylst/ml/__init__.py
+
+from .ml import train_regression_model
+from .ml import train_classification_model
+from .ml import train_test_split
+from .ml import accuracy_score
+
diff --git a/pylst/ml/ml.py b/pylst/ml/ml.py
new file mode 100644
index 0000000..a69d352
--- /dev/null
+++ b/pylst/ml/ml.py
@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Azad Rasul
+"""
+from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
+from sklearn.model_selection import train_test_split
+from sklearn.metrics import accuracy_score
+
+def train_regression_model(X_train, y_train, n_estimators=100):
+    """
+    Trains a Random Forest Regressor model.
+
+    Parameters:
+    - X_train: Training features
+    - y_train: Training labels
+    - n_estimators: Number of trees in the forest (default is 100)
+
+    Returns:
+    - Trained Random Forest Regressor model
+    """
+    # Initialize the Random Forest Regressor with the specified number of trees
+    regressor = RandomForestRegressor(n_estimators=n_estimators)
+    
+    # Train the model using the training data
+    regressor.fit(X_train, y_train)
+    
+    return regressor  # Return the trained model
+
+def train_classification_model(X_train, y_train, n_estimators=100):
+    """
+    Trains a Random Forest Classifier model.
+
+    Parameters:
+    - X_train: Training features
+    - y_train: Training labels
+    - n_estimators: Number of trees in the forest (default is 100)
+
+    Returns:
+    - Trained Random Forest Classifier model
+    """
+    # Initialize the Random Forest Classifier with the specified number of trees
+    classifier = RandomForestClassifier(n_estimators=n_estimators)
+    
+    # Train the classifier using the training data
+    classifier.fit(X_train, y_train)
+    
+    return classifier  # Return the trained model
+
+# Importing functions from sklearn for further use
+# train_test_split: Function to split data into training and testing sets
+# predict: Function to make predictions using RandomForestClassifier
+# accuracy_score: Function to calculate the accuracy of classification models
+train_test_split = train_test_split
+predict = RandomForestClassifier().predict
+accuracy_score = accuracy_score
diff --git a/pylst/normalization/__init__.py b/pylst/normalization/__init__.py
new file mode 100644
index 0000000..6dd953f
--- /dev/null
+++ b/pylst/normalization/__init__.py
@@ -0,0 +1,3 @@
+# pylst/normalization/__init__.py
+
+from .normalizer import nrs
diff --git a/pylst/normalization/normalizer.py b/pylst/normalization/normalizer.py
new file mode 100644
index 0000000..fac5cfe
--- /dev/null
+++ b/pylst/normalization/normalizer.py
@@ -0,0 +1,22 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Jan  7 14:08:35 2024
+
+@author: Azad Rasul
+"""
+
+# pylst/normalization/normalizer.py
+import numpy as np
+
+def nrs(lst):
+    # Compute NR:
+    NR = lst / np.sqrt(np.sum(np.square(lst)))
+    
+    # Multiply NR by the ratio of means
+    NRS = NR * (np.mean(lst) / np.mean(NR))
+    
+    return NRS
+    
+        
+        
+    
\ No newline at end of file
diff --git a/pylst/open/__init__.py b/pylst/open/__init__.py
new file mode 100644
index 0000000..dff86b8
--- /dev/null
+++ b/pylst/open/__init__.py
@@ -0,0 +1,6 @@
+# pylst/open/__init__.py
+
+from .opener import Opener
+from .remopener import RemoteOpener
+from .landsat_downloader import download_images
+
diff --git a/pylst/open/landsat_downloader.py b/pylst/open/landsat_downloader.py
new file mode 100644
index 0000000..4be3cef
--- /dev/null
+++ b/pylst/open/landsat_downloader.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Azad Rasul
+"""
+# pylst/open/landsat_downloader.py
+
+import ee
+from datetime import datetime
+
+def download_images(start_date, end_date, max_cloud_cover=20, region=None):
+    """
+    Download Landsat images from Google Earth Engine and export to Google Drive.
+
+    Parameters:
+    - start_date: Start date (formatted as 'YYYY-MM-DD').
+    - end_date: End date (formatted as 'YYYY-MM-DD').
+    - max_cloud_cover: Maximum allowed cloud cover percentage.
+    - region: A GeoJSON-like geometry specifying the region of interest (optional).
+
+    Returns:
+    - None
+    """
+    # Initialize the Earth Engine API
+    ee.Initialize()
+
+    # Define a time range
+    start_date = datetime.strptime(start_date, '%Y-%m-%d')
+    end_date = datetime.strptime(end_date, '%Y-%m-%d')
+
+    # Filter Landsat collection based on specified criteria
+    landsat_collection = (ee.ImageCollection('LANDSAT/LC08/C01/T1')
+                         .filterDate(start_date, end_date)
+                         .filterMetadata('CLOUD_COVER', 'less_than', max_cloud_cover)
+                         )
+
+    if region:
+        # Filter by region of interest
+        landsat_collection = landsat_collection.filterBounds(ee.Geometry.Polygon(region))
+
+    # Print the number of images in the collection
+    print("Number of images in the collection:", landsat_collection.size().getInfo())
+
+    # Iterate over the images and export to Google Drive
+    images = landsat_collection.toList(landsat_collection.size())
+    for i in range(images.size().getInfo()):
+        image = ee.Image(images.get(i))
+
+        # Get the Landsat ID
+        landsat_id = image.get('LANDSAT_ID').getInfo()
+
+        # Define export parameters
+        export_params = {
+            'image': image,
+            'description': "",  # Empty description
+            'folder': 'GEE_exports',  # Specify the Google Drive folder
+            'fileNamePrefix': f"Landsat_image_{i+1}",  # Use a prefix and system index
+            'scale': 30,
+            }            
+
+
+        # If a region is specified, add it to export parameters
+        if region:
+            export_params['region'] = region
+
+        # Export the image to Google Drive
+        task = ee.batch.Export.image.toDrive(**export_params)
+        task.start()
+        print(f"Export task {i+1} started:", task.id)
diff --git a/pylst/open/opener.py b/pylst/open/opener.py
new file mode 100644
index 0000000..5232044
--- /dev/null
+++ b/pylst/open/opener.py
@@ -0,0 +1,22 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Azad Rasul
+"""
+class Opener:
+    @staticmethod
+    def open(path, *args, **kwargs):
+        # Simplified logic to open image
+        try:
+            file = open(path, 'rb')
+            data = file.read()  # Read the file content
+            print(f"Opened file: {path}")
+            return file  # Return the file object
+        except FileNotFoundError:
+            print(f"File not found: {path}")
+            return None  # Return None if the file is not found
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, exc_type, exc_value, traceback):
+        pass
diff --git a/pylst/open/remopener.py b/pylst/open/remopener.py
new file mode 100644
index 0000000..a6aa906
--- /dev/null
+++ b/pylst/open/remopener.py
@@ -0,0 +1,22 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Azad Rasul
+"""
+# pylst/open/remopener.py
+import rasterio
+from io import BytesIO
+
+class RemoteOpener:
+    @staticmethod
+    def open(remote_url, *args, **kwargs):
+        try:
+            # Open the remote raster file
+            with rasterio.Env():
+                with rasterio.open(remote_url) as src:
+                    # Read the content of the raster file
+                    content = src.read()
+                    print(f"Successfully opened remote file: {remote_url}")
+                    return BytesIO(content.tobytes())
+        except rasterio.errors.RasterioError as e:
+            print(f"Failed to open remote file: {remote_url}. Error: {e}")
+            return None
diff --git a/pylst/pylst.py b/pylst/pylst.py
index dd0b80e..deddeb7 100644
--- a/pylst/pylst.py
+++ b/pylst/pylst.py
@@ -1 +1,235 @@
-"""Main module."""
+import numpy as np
+
+from .temperature import default_algorithms as temperature_algorithms
+from .emissivity import default_algorithms as emissivity_algorithms
+from .temperature import BrightnessTemperatureLandsat
+from .runner import Runner
+from .utils import compute_ndvi
+from .exceptions import *
+
+# Credit to the original author: Oladimeji Mudele (pylandtemp)
+# Original source: https://github.com/pylandtemp/pylandtemp
+
+def split_window(
+    landsat_band_10: np.ndarray,
+    landsat_band_11: np.ndarray,
+    landsat_band_4: np.ndarray,
+    landsat_band_5: np.ndarray,
+    lst_method: str,
+    emissivity_method: str,
+    unit: str = "kelvin",
+) -> np.ndarray:
+    """Provides an interface to compute land surface temperature
+        from landsat 8 imagery using split window method
+
+    Args:
+        landsat_band_10 (np.ndarray): Band 10 of the Landsat 8 image
+        landsat_band_11 (np.ndarray): Band 11 of the landsat 8 image
+        landsat_band_4 (np.ndarray): Band 4 of the landsat 8 image (Red band)
+        landsat_band_5 (np.ndarray): Band 5 of the landsat 8 image (Near-Infrared band)
+
+        lst_method (str): provide one of the valid split window method for computing land surface temperature
+                        Valid methods to add include:
+                        'jiminez-munoz': Jiminez-Munoz et al, 2008
+                        'kerr': Kerr Y et al, 2004
+                        'mc-millin': McMillin L. M. , 1975
+                        'price': Price J. C., 1984
+                        'sobrino-1993': Sobrino J. A. et al, 1993
+                        'coll-caselles': Coll C. et al, 1997
+
+        emissivity_method (str): provide one of the valid split window method for computing land surface emissivity.
+                                         Valid methods to add include:
+                                        'advan': Avdan Ugur et al, 2016
+                                        'xiaolei':  Xiaolei Yu et al, 2014
+                                        'gopinadh': Gopinadh Rongali et al 2018
+
+        unit (str, optional): 'kelvin' or 'celsius'. Defaults to 'kelvin'.
+
+    Returns:
+        np.ndarray: Land surface temperature (numpy array)
+    """
+
+    if not (
+        landsat_band_10.shape
+        == landsat_band_11.shape
+        == landsat_band_5.shape
+        == landsat_band_4.shape
+    ):
+        raise InputShapesNotEqual(
+            f"Shapes of input images should be equal: {landsat_band_10.shape}, {landsat_band_5.shape}, {landsat_band_4.shape}"
+        )
+
+    mask = landsat_band_10 == 0
+    ndvi_image = ndvi(landsat_band_5, landsat_band_4, mask)
+
+    brightness_temp_10, brightness_temp_11 = brightness_temperature(
+        landsat_band_10, landsat_band_11=landsat_band_11, mask=mask
+    )
+
+    emissivity_10, emissivity_11 = Runner(algorithms=emissivity_algorithms)(
+        emissivity_method, ndvi=ndvi_image, red_band=landsat_band_4
+    )
+
+    lst_image = Runner(algorithms=temperature_algorithms.split_window)(
+        lst_method,
+        emissivity_10=emissivity_10,
+        emissivity_11=emissivity_11,
+        brightness_temperature_10=brightness_temp_10,
+        brightness_temperature_11=brightness_temp_11,
+        mask=mask,
+        ndvi=ndvi_image,
+    )
+    return lst_image
+
+
+def single_window(
+    landsat_band_10: np.ndarray,
+    landsat_band_4: np.ndarray,
+    landsat_band_5: np.ndarray,
+    lst_method: str = "mono-window",
+    emissivity_method: str = "avdan",
+    unit: str = "kelvin",
+) -> np.ndarray:
+    """Provides an interface to compute land surface temperature
+        from landsat 8 imagery using single window method
+
+    Args:
+        landsat_band_10 (np.ndarray): Band 10 of the Landsat 8 image
+        landsat_band_4 (np.ndarray): Band 4 of the Landsat 8 image (Red band)
+        landsat_band_5 (np.ndarray): Band 4 of the Landsat 8 image (Near-Infrared band)
+
+        lst_method (str, optional): Defaults to 'mono-window'.
+                                    Valid methods to add include:
+                                    'mono-window': Avdan Ugur et al, 2016
+
+        emissivity_method (str, optional): provide one of the valid split window method for computing
+                                        land surface emissivity.
+
+                                        Defaults to 'avdan'.
+                                         Valid methods to add include:
+                                        'advan': Avdan Ugur et al, 2016
+                                        'xiaolei':  Xiaolei Yu et al, 2014
+                                        'gopinadh': Gopinadh Rongali et al 2018
+
+        unit (str, optional): 'celsius' or 'kelvin'. Defaults to 'kelvin'.
+
+    Returns:
+        np.ndarray: Land surface temperature (numpy array)
+    """
+
+    if not landsat_band_10.shape == landsat_band_5.shape == landsat_band_4.shape:
+        raise InputShapesNotEqual(
+            f"Shapes of input images should be equal: {landsat_band_10.shape}, {landsat_band_5.shape}, {landsat_band_4.shape}"
+        )
+
+    mask = landsat_band_10 == 0
+    ndvi_image = ndvi(landsat_band_5, landsat_band_4, mask)
+    brightness_temp_10, _ = brightness_temperature(landsat_band_10, mask=mask)
+
+    emissivity_10, _ = Runner(algorithms=emissivity_algorithms)(
+        emissivity_method, ndvi=ndvi_image, red_band=landsat_band_4
+    )
+
+    lst_image = Runner(algorithms=temperature_algorithms.single_window)(
+        lst_method,
+        emissivity_10=emissivity_10,
+        brightness_temperature_10=brightness_temp_10,
+        mask=mask,
+        ndvi=ndvi_image,
+    )
+    return lst_image
+
+
+def emissivity(
+    ndvi_image: np.ndarray,
+    landsat_band_4: np.ndarray = None,
+    emissivity_method: str = "avdan",
+):
+    """Provides an interface to compute land surface emissivity
+        from landsat 8 imagery
+
+    Args:
+        ndvi_image (np.ndarray): Normalized difference vegetation index image
+
+        landsat_band_4 (None or np.ndarray, optional): red band image. Defaults to None.
+                                                        Can be None except when emissivity_method = 'xiaolei'
+
+        emissivity_method (str, optional): provide one of the valid split window method for computing land surface emissivity.
+                                            Defaults to 'avdan'.
+                                         Valid methods to add include:
+                                        'advan': Avdan Ugur et al, 2016
+                                        'xiaolei':  Xiaolei Yu et al, 2014
+
+    Returns:
+        np.ndarray: Emissivity numpy array
+    """
+    if not ndvi_image.shape == landsat_band_4.shape:
+        raise InputShapesNotEqual(
+            f"Shapes of input images should be equal: {ndvi_image.shape}, {landsat_band_4.shape}"
+        )
+    if emissivity_method == "xiaolei" and landsat_band_4 is None:
+        raise ValueError(
+            f"The red band (landsat_band_4) has to be provided if {emissivity_method} is to be used"
+        )
+
+    emissivity_10, emissivity_11 = Runner(algorithms=emissivity_algorithms)(
+        emissivity_method, ndvi=ndvi_image, red_band=landsat_band_4
+    )
+    return emissivity_10, emissivity_11
+
+
+def ndvi(landsat_band_5: np.ndarray, landsat_band_4: np.ndarray, mask: np.ndarray):
+    """Computes the NDVI given bands 4 and 5 of Landsat 8 image
+
+    Args:
+        landsat_band_5 (np.ndarray): Band 5 of landsat 8 image
+        landsat_band_4 (np.ndarray): Band 4 of landsat 8 image
+        mask (np.ndarray[bool]): output is NaN where Mask == True
+
+    Returns:
+        np.ndarray: NVDI numpy array
+    """
+    if not landsat_band_5.shape == landsat_band_4.shape:
+        raise InputShapesNotEqual(
+            f"Shapes of input images should be equal: {landsat_band_5.shape}, {landsat_band_4.shape}"
+        )
+
+    if mask.dtype != "bool":
+        raise InvalidMaskError(
+            f"Image passed in as 'mask' must be a numpy array with bool dtype values"
+        )
+    return compute_ndvi(landsat_band_5, landsat_band_4, mask=mask)
+
+
+def brightness_temperature(
+    landsat_band_10: np.ndarray,
+    landsat_band_11: np.ndarray = None,
+    mask: np.ndarray = None,
+):
+    """Compute brightness temperature
+
+    Args:
+        landsat_band_10 (np.ndarray): Band 10 of landsat 8 image
+        landsat_band_11 (np.ndarray): Band 11 of landsat 8 image. Defaults to None.
+        mask (np.ndarray[bool]): output is NaN where Mask == True. Defaults to None.
+
+    Returns:
+        np.ndarray: Brightness temperature numpy array
+    """
+    if (
+        landsat_band_11 is not None
+        and not landsat_band_10.shape == landsat_band_11.shape
+    ):
+        raise InputShapesNotEqual(
+            f"Shapes of input images should be equal: {landsat_band_10.shape}, {landsat_band_11.shape}"
+        )
+
+    if mask.dtype != "bool":
+        raise InvalidMaskError(
+            f"image passed in as 'mask' must be a numpy array with bool dtype values"
+        )
+
+    brightness_temp_10, brightness_temp_11 = BrightnessTemperatureLandsat()(
+        landsat_band_10, landsat_band_11, mask=mask
+    )
+    return brightness_temp_10, brightness_temp_11
diff --git a/pylst/runner.py b/pylst/runner.py
new file mode 100644
index 0000000..8facc78
--- /dev/null
+++ b/pylst/runner.py
@@ -0,0 +1,20 @@
+class Runner:
+    def __init__(self, algorithms):
+        """
+        Task and method agnostic runner
+
+        args:
+        algorithms (dict): dictionary that maps defined keys to concrete implementations of the called algorithm.
+        """
+        self.algorithms = algorithms
+
+    def __call__(self, method, **kwargs):
+        compute_algorithm = self._get_algorithm(method)
+        return compute_algorithm()(**kwargs)
+
+    def _get_algorithm(self, algo):
+        if algo not in self.algorithms:
+            raise ValueError(
+                f"Requested method not implemented. Choose among available methods: {list(self.algorithms.values())}"
+            )
+        return self.algorithms[algo]
diff --git a/pylst/spatial_analysis/__init__.py b/pylst/spatial_analysis/__init__.py
new file mode 100644
index 0000000..a2eccfe
--- /dev/null
+++ b/pylst/spatial_analysis/__init__.py
@@ -0,0 +1,6 @@
+# __init__.py
+
+# Importing the calculate_zs function from pylst.spatial_analysis.zonstat
+# This function is used for spatial analysis and zone statistics.
+from pylst.spatial_analysis.zonstat import calculate_zs
+from pylst.spatial_analysis.changedet import analyze_images
\ No newline at end of file
diff --git a/pylst/spatial_analysis/changedet.py b/pylst/spatial_analysis/changedet.py
new file mode 100644
index 0000000..61dcb3a
--- /dev/null
+++ b/pylst/spatial_analysis/changedet.py
@@ -0,0 +1,59 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Azad Rasul
+"""
+# Import necessary libraries for image analysis
+import rasterio
+import numpy as np
+
+def analyze_images(image_path1, image_path2, output_path):
+    """
+    Analyzes two images and saves the difference as a GeoTIFF.
+
+    Parameters:
+    - image_path1 (str): File path to the first image.
+    - image_path2 (str): File path to the second image.
+    - output_path (str): File path to save the analysis result as GeoTIFF.
+
+    Returns:
+    - result (dict): A dictionary containing the analysis results or None if an error occurs.
+    """
+    try:
+        # Open the images using rasterio
+        with rasterio.open(image_path1) as src1, rasterio.open(image_path2) as src2:
+            # Read the image data
+            image1 = src1.read(1)  # Assuming a single band for simplicity
+            image2 = src2.read(1)  # Assuming a single band for simplicity
+
+            # Check if the shapes of the two images are compatible
+            if image1.shape != image2.shape:
+                raise ValueError("Images have different shapes and cannot be subtracted.")
+
+            # Perform any desired analysis on the images
+            # For example, you can calculate the difference between the two images
+            image_difference = np.subtract(image2, image1)
+
+            # You can add more analysis steps here based on your specific requirements
+
+            # Store the results in a dictionary
+            result = {
+                "image_difference": image_difference,
+                # Add more result keys and values as needed
+            }
+
+            # Get metadata from one of the input images (assuming both have the same properties)
+            meta = src1.meta.copy()
+
+            # Update metadata for the output GeoTIFF
+            meta.update(driver='GTiff', count=1, dtype='float32')
+
+            # Create a new GeoTIFF file for the image_difference
+            with rasterio.open(output_path, 'w', **meta) as dst:
+                dst.write(image_difference.astype('float32'), 1)
+
+            return result
+
+    except Exception as e:
+        # Handle exceptions such as file not found or invalid file format
+        print(f"Error during image analysis: {e}")
+        return None
diff --git a/pylst/spatial_analysis/zonstat.py b/pylst/spatial_analysis/zonstat.py
new file mode 100644
index 0000000..9bf94e1
--- /dev/null
+++ b/pylst/spatial_analysis/zonstat.py
@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Azad Rasul
+"""
+# Import necessary libraries
+import pandas as pd
+from rasterstats import zonal_stats
+
+def calculate_zs(shapefile_path, raster_path):
+    """
+    Calculate zonal statistics for a given shapefile and raster.
+
+    Parameters:
+    - shapefile_path (str): Path to the shapefile containing the zones of interest.
+    - raster_path (str): Path to the raster file for which zonal statistics are calculated.
+
+    Returns:
+    - df (pd.DataFrame): Pandas DataFrame containing the calculated zonal statistics.
+    """
+
+    # Extract zonal statistics using the rasterstats library
+    stats = zonal_stats(shapefile_path, raster_path, stats=["min", "max", "mean", "median", "std", "count", "range", "sum",
+                                                            "median", "majority", "minority", "unique", "nodata"], geojson_out=False)
+    
+    # Convert the results to a pandas DataFrame for easier analysis and manipulation
+    df = pd.DataFrame(stats)
+
+    # Save the DataFrame as a CSV file for further reference or sharing
+    df.to_csv("zonal_stats.csv", index=True)
+
+    # Return the DataFrame containing the zonal statistics
+    return df
diff --git a/pylst/temperature/__init__.py b/pylst/temperature/__init__.py
new file mode 100644
index 0000000..2d9ead6
--- /dev/null
+++ b/pylst/temperature/__init__.py
@@ -0,0 +1,19 @@
+# Importing necessary components from different modules
+
+# Importing default algorithms for temperature computation
+from .temperature import default_algorithms
+
+# Importing the BrightnessTemperatureLandsat class for brightness temperature calculation
+from .brightness_temperature import BrightnessTemperatureLandsat
+
+# Importing the MonoWindowLST algorithm for mono-window land surface temperature calculation
+from .algorithms.mono_window import MonoWindowLST
+
+# Importing split window algorithms for land surface temperature calculation
+from .algorithms.split_window.algorithms import (
+    SplitWindowJiminezMunozLST,
+    SplitWindowKerrLST,
+    SplitWindowMcMillinLST,
+    SplitWindowPriceLST,
+    SplitWindowSobrino1993LST,
+)
\ No newline at end of file
diff --git a/pylst/temperature/algorithms/__init__.py b/pylst/temperature/algorithms/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/pylst/temperature/algorithms/mono_window.py b/pylst/temperature/algorithms/mono_window.py
new file mode 100644
index 0000000..b144d7c
--- /dev/null
+++ b/pylst/temperature/algorithms/mono_window.py
@@ -0,0 +1,75 @@
+import numpy as np
+from pylst.exceptions import assert_required_keywords_provided
+
+# Credit to the original author: Oladimeji Mudele (pylandtemp)
+# Original source: https://github.com/pylandtemp/pylandtemp
+
+class MonoWindowLST:
+    def __init__(self):
+        """
+        Method reference:
+        Avdan, Ugur, and Gordana Jovanovska. "Algorithm for automated mapping of land surface
+        temperature using LANDSAT 8 satellite data." Journal of sensors 2016 (2016).
+
+        Initializes the Mono-window Land Surface Temperature (LST) computation algorithm.
+
+        Attributes:
+            max_earth_temp (float): Maximum possible Earth temperature in Kelvin.
+        """
+        self.max_earth_temp = 273.15 + 56.7
+
+    def __call__(self, **kwargs) -> np.ndarray:
+        """
+        Compute land surface temperature using the Mono-window algorithm.
+
+        Args:
+            emissivity_10 (np.ndarray): Emissivity image obtained for band 10
+            brightness_temperature_10 (np.ndarray): Brightness temperature image obtained for band 10
+            mask (np.ndarray[bool]): Mask image. Output will have NaN value where mask is True.
+
+        Returns:
+            np.ndarray: Land surface temperature image
+        """
+        self._validate_inputs(**kwargs)
+        return self._compute_lst_mono_window(**kwargs)
+
+    def _validate_inputs(self, **kwargs):
+        """
+        Validate input parameters for the Mono-window LST computation.
+
+        Args:
+            **kwargs: Keyword arguments containing necessary input images.
+
+        Raises:
+            ValueError: If input images have different sizes.
+        """
+        
+        required_keywords = ["brightness_temperature_10", "emissivity_10", "mask"]
+        assert_required_keywords_provided(required_keywords, **kwargs)
+        temperature_band = kwargs["brightness_temperature_10"]
+        emissivity = kwargs["emissivity_10"]
+        mask = kwargs["mask"]
+
+        if not (mask.shape == temperature_band.shape == emissivity.shape):
+            raise ValueError("Input images must be of the same size/shape")
+
+    def _compute_lst_mono_window(self, **kwargs) -> np.ndarray:
+        """
+        Compute land surface temperature using the Mono-window algorithm.
+
+        Args:
+            **kwargs: Keyword arguments containing necessary input images.
+
+        Returns:
+            np.ndarray: Land surface temperature image.
+        """
+        temperature_band = kwargs["brightness_temperature_10"]
+        emissivity = kwargs["emissivity_10"]
+        mask = kwargs["mask"]
+
+        land_surface_temp = temperature_band / (
+            1 + (((0.0000115 * temperature_band) / 14380) * np.log(emissivity))
+        )
+        land_surface_temp[mask] = np.nan
+        land_surface_temp[land_surface_temp > self.max_earth_temp] = np.nan
+        return land_surface_temp
diff --git a/pylst/temperature/algorithms/split_window/__init__.py b/pylst/temperature/algorithms/split_window/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/pylst/temperature/algorithms/split_window/__pycache__/__init__.cpython-39.pyc b/pylst/temperature/algorithms/split_window/__pycache__/__init__.cpython-39.pyc
new file mode 100644
index 0000000..3c8d8a7
Binary files /dev/null and b/pylst/temperature/algorithms/split_window/__pycache__/__init__.cpython-39.pyc differ
diff --git a/pylst/temperature/algorithms/split_window/__pycache__/algorithms.cpython-39.pyc b/pylst/temperature/algorithms/split_window/__pycache__/algorithms.cpython-39.pyc
new file mode 100644
index 0000000..b692ed9
Binary files /dev/null and b/pylst/temperature/algorithms/split_window/__pycache__/algorithms.cpython-39.pyc differ
diff --git a/pylst/temperature/algorithms/split_window/algorithms.py b/pylst/temperature/algorithms/split_window/algorithms.py
new file mode 100644
index 0000000..c47b1e6
--- /dev/null
+++ b/pylst/temperature/algorithms/split_window/algorithms.py
@@ -0,0 +1,254 @@
+import numpy as np
+from pylst.exceptions import assert_required_keywords_provided
+from pylst.utils import fractional_vegetation_cover
+
+# Credit to the original author: Oladimeji Mudele (pylandtemp)
+# Original source: https://github.com/pylandtemp/pylandtemp
+
+class SplitWindowParentLST:
+    def __init__(self):
+        # Maximum Earth temperature for clipping the calculated LST
+        self.max_earth_temp = 273.15 + 56.7
+
+    def __call__(self, **kwargs) -> np.ndarray:
+        # Calculate LST and clip values exceeding the maximum Earth temperature
+        lst = self._compute_lst(**kwargs)
+        lst[lst > self.max_earth_temp] = np.nan
+        return lst
+
+    def _compute_lst(self, **kwargs):
+        # Abstract method for computing LST, to be implemented in derived classes
+        raise NotImplementedError("Concrete method yet to be implemented")
+
+
+class SplitWindowJiminezMunozLST(SplitWindowParentLST):
+    cwv = 0.013
+
+    def _compute_lst(self, **kwargs) -> np.ndarray:
+        # Concrete implementation for Jiminez-Munoz method
+        required_keywords = [
+            "emissivity_10",
+            "emissivity_11",
+            "brightness_temperature_10",
+            "brightness_temperature_11",
+            "mask",
+        ]
+        assert_required_keywords_provided(required_keywords, **kwargs)
+        
+        # Extracting required input parameters
+        tb_10 = kwargs["brightness_temperature_10"]
+        tb_11 = kwargs["brightness_temperature_11"]
+        emissivity_10 = kwargs["emissivity_10"]
+        emissivity_11 = kwargs["emissivity_11"]
+        mask = kwargs["mask"]
+        
+        # Computing LST using Jiminez-Munoz method
+        mean_e = (emissivity_10 + emissivity_11) / 2
+        diff_e = emissivity_10 - emissivity_11
+        diff_tb = tb_10 - tb_11
+
+        lst = (
+            tb_10
+            + (1.387 * diff_tb)
+            + (0.183 * (diff_tb ** 2))
+            - 0.268
+            + ((54.3 - (2.238 * self.cwv)) * (1 - mean_e))
+            + ((-129.2 + (16.4 * self.cwv)) * diff_e)
+        )
+        lst[mask] = np.nan
+        return lst
+
+
+class SplitWindowKerrLST(SplitWindowParentLST):
+    """
+    Method reference:
+
+    Kerr Y, Lagouarde J, Nerry F, Ottlé C (2004) Land surface temperature
+    retrieval techniques and applications: case of the AVHRR. Thermal remote
+    sensing in land surface processing. CRC Press, New York, pp 55–131
+    """
+
+    def _compute_lst(self, **kwargs) -> np.ndarray:
+        """
+        Computes the Land Surface Temperature (LST) using the Split Window method based on Kerr et al. (2004).
+
+        Args:
+            brightness_temperature_10 (np.ndarray): Brightness temperature image obtained for band 10
+            brightness_temperature_11 (np.ndarray): Brightness temperature image obtained for band 11
+            ndvi (np.ndarray): NDVI image
+            mask (np.ndarray[bool]): Mask image. Output will have NaN value where the mask is True.
+
+        Returns:
+            np.ndarray: Land surface temperature image
+        """
+        # Ensure that required keyword arguments are provided
+        required_keywords = [
+            "brightness_temperature_10",
+            "brightness_temperature_11",
+            "ndvi",
+            "mask",
+        ]
+        assert_required_keywords_provided(required_keywords, **kwargs)
+        
+        # Retrieve required input data
+        tb_10 = kwargs["brightness_temperature_10"]
+        tb_11 = kwargs["brightness_temperature_11"]
+        ndvi = kwargs["ndvi"]
+        mask = kwargs["mask"]
+        
+        # Calculate fractional vegetation cover (pv) using NDVI
+        pv = fractional_vegetation_cover(ndvi)
+
+        # Compute Land Surface Temperature (LST) using the Kerr et al. (2004) formula
+        lst = (
+            (tb_10 * ((0.5 * pv) + 3.1))
+            + (tb_11 * ((-0.5 * pv) - 2.1))
+            - ((5.5 * pv) + 3.1)
+        )
+        
+        # Mask out invalid values based on the provided mask
+        lst[mask] = np.nan
+        
+        return lst
+
+
+class SplitWindowMcMillinLST(SplitWindowParentLST):
+    """
+    Method reference:
+    McMillin LM (1975) Estimation of sea surface temperatures
+    from two infrared window measurements with different absorption.
+    J Geophys Res 80(36):5113–5117. https://doi.org/10.1029/JC080i036p05113
+    """
+
+    def _compute_lst(self, **kwargs) -> np.ndarray:
+        """
+        Computes the land surface temperature using McMillin's method.
+
+        Args:
+            **brightness_temperature_10 (np.ndarray): Brightness temperature image obtained for band 10
+            **brightness_temperature_11 (np.ndarray): Brightness temperature image obtained for band 11
+            **mask (np.ndarray[bool]): Mask image. Output will have NaN values where the mask is True.
+
+        Returns:
+            np.ndarray: Land surface temperature image
+        """
+        # Check if required keywords are provided
+        required_keywords = [
+            "brightness_temperature_10",
+            "brightness_temperature_11",
+            "mask",
+        ]
+        assert_required_keywords_provided(required_keywords, **kwargs)
+
+        # Extracting input variables
+        tb_10 = kwargs["brightness_temperature_10"]
+        tb_11 = kwargs["brightness_temperature_11"]
+        mask = kwargs["mask"]
+
+        # McMillin's land surface temperature computation
+        lst = (1.035 * tb_10) + (3.046 * (tb_10 - tb_11)) - 10.93
+
+        # Masking areas where the mask is True
+        lst[mask] = np.nan
+
+        return lst
+
+
+
+class SplitWindowPriceLST(SplitWindowParentLST):
+    """
+    Method reference:
+    Price JC (1984) Land surface temperature measurements from the split window
+    channels of the NOAA advanced very high-resolution radiometer. J Geophys Res 89:7231–7237.
+    https://doi.org/10.1029/JD089iD05p07231
+    """
+
+    def _compute_lst(self, **kwargs) -> np.ndarray:
+        """
+        Computes Land Surface Temperature (LST) using the split window Price method.
+
+        Args:
+            **kwargs (dict): Keyword arguments containing necessary input arrays.
+                             Required keywords: "emissivity_10", "emissivity_11",
+                             "brightness_temperature_10", "brightness_temperature_11", "mask".
+
+        Returns:
+            np.ndarray: Land Surface Temperature (LST) array.
+        """
+        required_keywords = [
+            "emissivity_10",
+            "emissivity_11",
+            "brightness_temperature_10",
+            "brightness_temperature_11",
+            "mask",
+        ]
+        assert_required_keywords_provided(required_keywords, **kwargs)
+
+        # Extracting required arrays from kwargs
+        tb_10 = kwargs["brightness_temperature_10"]
+        tb_11 = kwargs["brightness_temperature_11"]
+        emm_10 = kwargs["emissivity_10"]
+        emm_11 = kwargs["emissivity_11"]
+        mask = kwargs["mask"]
+
+        # Computing LST using the Price method formula
+        lst = (tb_10 + 3.33 * (tb_10 - tb_11)) * ((5.5 - emm_10) / 4.5) + (
+            0.75 * tb_11 * (emm_10 - emm_11)
+        )
+        
+        # Masking invalid values in the result array
+        lst[mask] = np.nan
+        return lst
+
+
+
+class SplitWindowSobrino1993LST(SplitWindowParentLST):
+    """
+    Method reference:
+
+    Sobrino JA, Caselles V, Coll C (1993) Theoretical split window algorithms
+    for determining the actual surface temperature. I1Nuovo Cimento 16:219–236.
+    https://doi.org/10.1007/BF02524225
+
+    This class implements the Split Window algorithm based on the Sobrino 1993 method.
+    The algorithm computes the land surface temperature (LST) using brightness temperatures
+    and emissivity values for Landsat 8 bands 10 and 11.
+    """
+
+    def _compute_lst(self, **kwargs) -> np.ndarray:
+        # Define the required keywords for the algorithm
+        required_keywords = [
+            "emissivity_10",
+            "emissivity_11",
+            "brightness_temperature_10",
+            "brightness_temperature_11",
+            "mask",
+        ]
+        # Ensure that all required keywords are provided
+        assert_required_keywords_provided(required_keywords, **kwargs)
+
+        # Retrieve necessary data from the provided keyword arguments
+        tb_10 = kwargs["brightness_temperature_10"]
+        tb_11 = kwargs["brightness_temperature_11"]
+        emm_10 = kwargs["emissivity_10"]
+        emm_11 = kwargs["emissivity_11"]
+        mask = kwargs["mask"]
+
+        # Calculate the temperature difference between tb_10 and tb_11
+        diff_tb = tb_10 - tb_11
+        # Calculate the emissivity difference between emm_10 and emm_11
+        diff_e = emm_10 - emm_11
+
+        # Compute the land surface temperature (LST) using the Sobrino 1993 formula
+        lst = (
+            tb_10
+            + (1.06 * diff_tb)
+            + (0.46 * diff_tb ** 2)
+            + (53 * (1 - emm_10))
+            - (53 * diff_e)
+        )
+        # Mask invalid pixels based on the provided mask, setting their values to NaN
+        lst[mask] = np.nan
+        return lst
+
+
diff --git a/pylst/temperature/brightness_temperature.py b/pylst/temperature/brightness_temperature.py
new file mode 100644
index 0000000..0bc875b
--- /dev/null
+++ b/pylst/temperature/brightness_temperature.py
@@ -0,0 +1,56 @@
+import numpy as np
+from .utils import compute_brightness_temperature
+
+# Credit to the original author: Oladimeji Mudele (pylandtemp)
+# Original source: https://github.com/pylandtemp/pylandtemp
+
+class BrightnessTemperatureLandsat:
+    def __init__(self):
+        # Constants for brightness temperature calculation
+        self.constants = {
+            "band_10": {"k1": 774.89, "k2": 1321.08},
+            "band_11": {"k1": 480.89, "k2": 1201.14},
+            "mult_factor": 0.0003342,
+            "add_factor": 0.1,
+        }
+
+    def __call__(self, band_10: np.ndarray, band_11: np.ndarray = None, mask=None) -> np.ndarray:
+        """
+        Calculate brightness temperature for Landsat bands 10 and 11.
+
+        Args:
+            band_10 (np.ndarray): Level 1 quantized and calibrated scaled Digital Numbers (DN) TIR band data for Band 10 landsat 8 data.
+            band_11 (np.ndarray): Level 1 quantized and calibrated scaled Digital Numbers (DN) TIR band data for Band 11 landsat 8 data.
+            mask (bool): Mask zero or NaN values. Defaults to True.
+
+        Returns:
+            Tuple(np.ndarray, np.ndarray): Band 10 brightness temperature, Band 11 brightness temperature.
+        """
+        # Compute brightness temperature for Band 10
+        tb_band_10 = self._compute_brightness_temp(band_10, "band_10", mask)
+        tb_band_11 = None
+
+        if band_11 is not None:
+            # Compute brightness temperature for band 11 if available
+            tb_band_11 = self._compute_brightness_temp(band_11, "band_11", mask)
+
+        return tb_band_10, tb_band_11
+
+    def _compute_brightness_temp(self, image: np.ndarray, band: str, mask: np.ndarray) -> np.ndarray:
+        """Converts image raw digital numbers to brightness temperature.
+
+        Args:
+            image (np.ndarray): Band image data.
+            band (str): Band identifier.
+            mask (n.ndarray[bool]): True for pixels to mask out.
+
+        Returns:
+            np.ndarray: Brightness temperature corrected image.
+        """
+        k1 = self.constants[band]["k1"]
+        k2 = self.constants[band]["k2"]
+        mult_factor = self.constants["mult_factor"]
+        add_factor = self.constants["add_factor"]
+        
+        # Use utility function to compute brightness temperature
+        return compute_brightness_temperature(image, mult_factor, add_factor, k1, k2, mask)
diff --git a/pylst/temperature/temperature.py b/pylst/temperature/temperature.py
new file mode 100644
index 0000000..6de734a
--- /dev/null
+++ b/pylst/temperature/temperature.py
@@ -0,0 +1,29 @@
+from collections import namedtuple
+from .algorithms.mono_window import MonoWindowLST
+from .algorithms.split_window.algorithms import (
+    SplitWindowJiminezMunozLST,
+    SplitWindowKerrLST,
+    SplitWindowMcMillinLST,
+    SplitWindowPriceLST,
+    SplitWindowSobrino1993LST,
+)
+
+# Define namedtuples to organize available algorithms
+SingleWindowAlgorithms = {
+    "mono-window": MonoWindowLST # Algorithm for single window method
+}
+
+SplitWindowAlgorithms = {
+    "jiminez-munoz": SplitWindowJiminezMunozLST,
+    "kerr": SplitWindowKerrLST,
+    "mc-millin": SplitWindowMcMillinLST,
+    "price": SplitWindowPriceLST,
+    "sobrino-1993": SplitWindowSobrino1993LST
+}
+
+# Organize algorithms using a namedtuple for clarity
+Algorithms = namedtuple("Algorithms", ("single_window", "split_window"))
+default_algorithms = Algorithms(SingleWindowAlgorithms, SplitWindowAlgorithms)
+
+# 'default_algorithms' contains available algorithms for both single and split window methods
+# Access them using 'default_algorithms.single_window' and 'default_algorithms.split_window'
diff --git a/pylst/temperature/utils.py b/pylst/temperature/utils.py
new file mode 100644
index 0000000..a269222
--- /dev/null
+++ b/pylst/temperature/utils.py
@@ -0,0 +1,37 @@
+import numpy as np
+
+def compute_brightness_temperature(
+    image: np.ndarray, M: float, A: float, k1: float, k2: float, mask: np.ndarray = None
+) -> np.ndarray:
+    """Converts Landsat thermal band Digital Numbers (DN) to brightness temperature.
+
+    Args:
+        image (np.ndarray): Level 1 quantized and calibrated scaled Digital Numbers
+                            (DN) TIR band data (e.g Band 10 landsat 8 data).
+        M (float): Band-specific multiplicative rescaling factor from the image
+                   folder metadata (RADIANCE_MULT_BAND_x, where x is the band number).
+        A (float): Band-specific additive rescaling factor from the image
+                   folder metadata (RADIANCE_ADD_BAND_x, where x is the band number).
+        k1 (float): Band-specific thermal conversion constant from the image
+                    folder metadata (K1_CONSTANT_BAND_x, where x is the thermal band number).
+        k2 (float): Band-specific thermal conversion constant from the image
+                    folder metadata (K2_CONSTANT_BAND_x, where x is the thermal band number).
+        mask (np.ndarray, optional): Boolean mask to identify NaN or irregular values
+                                     for masking purposes.
+
+    Returns:
+        np.ndarray: Brightness temperature-corrected Landsat image.
+    """
+    
+    # Convert raw digital numbers to top-of-atmoshere radiance
+    toa_radiance = (M * image) + A
+    
+    #  Compute brightness temperature using the Plank radiation law
+    brightness_temp = k2 / np.log((k1 / toa_radiance) + 1)
+    
+    
+    # Apply optional mask to handle NaN or irregular values
+    if mask is not None:
+        brightness_temp[mask] = np.nan
+
+    return brightness_temp
diff --git a/pylst/utils.py b/pylst/utils.py
new file mode 100644
index 0000000..d416309
--- /dev/null
+++ b/pylst/utils.py
@@ -0,0 +1,92 @@
+import numpy as np
+
+# Credit to the original author: Oladimeji Mudele (pylandtemp)
+# Original source: https://github.com/pylandtemp/pylandtemp
+
+def generate_mask(image: np.ndarray) -> np.ndarray:
+    """
+    Return a bool array masking 0 and NaN values as False and others as True
+
+    Args:
+        image (np.ndarray): Single-band image which is True where we do not want to mask and False where we want to mask.
+    """
+    zero_mask = image != 0
+    nan_mask = image != np.nan
+    mask_true = np.logical_and(zero_mask, nan_mask)
+    return mask_true
+
+
+def compute_ndvi(
+    nir: np.ndarray, red: np.ndarray, eps: float = 1e-15, mask=None
+) -> np.ndarray:
+    """Takes the near infrared and red bands of an optical satellite image as input and returns the ndvi: normalized difference vegetation index
+
+    Args:
+        nir (np.ndarray): Near-infrared band image
+        red (np.ndarray): Red-band image
+        eps (float): Epsilon to avoid ZeroDivisionError in numpy
+        use_mask (bool): If True, mask NaN and 0 values in input images.
+
+    Returns:
+        np.ndarray: Normalized difference vegetation index
+    """
+    ndvi = (nir - red) / (nir + red + eps)
+    ndvi[abs(ndvi) > 1] = np.nan
+    if mask is not None:
+        ndvi[mask] = np.nan
+    return ndvi
+
+
+def fractional_vegetation_cover(ndvi: np.ndarray) -> np.ndarray:
+    """Computes the fractinal vegetation cover matrix
+
+    Args:
+        ndvi (np.ndarray):  Normalized difference vegetation index (m x n)
+    Returns:
+        np.ndarray: Fractional vegetation cover
+    """
+    if len(ndvi.shape) != 2:
+        raise ValueError("NDVI image should be 2-dimensional")
+    return ((ndvi - 0.2) / (0.5 - 0.2)) ** 2
+
+
+def cavity_effect(
+    emissivity_veg: float,
+    emissivity_soil: float,
+    fractional_vegetation_cover: np.ndarray,
+    geometrical_factor: float = 0.55,
+) -> np.ndarray:
+    """Compute the cavity effect matrix
+
+    Args:
+        emissivity_veg (float): value of vegetation emissivity
+        emissivity_soil (float): value of soil emissivity
+        fractional_vegetation_cover (np.ndarray): Fractional vegetation cover image
+        geometrical_factor (float, optional): Geometric factor. Defaults to 0.55.
+
+    Returns:
+        np.ndarray: Cavity effect numpy array
+    """
+    to_return = (
+        (1 - emissivity_soil)
+        * emissivity_veg
+        * geometrical_factor
+        * (1 - fractional_vegetation_cover)
+    )
+    return to_return
+
+
+def rescale_band(
+    image: np.ndarray, mult: float = 2e-05, add: float = 0.1
+) -> np.ndarray:
+    """rescales the image band
+
+    Args:
+        image (np.ndarray): Band 1 - 9, or non Thermal IR bands of the satellite image.
+        mult (float, optional): Multiplicative factor. Defaults to 2e-05.
+        add (float, optional): Additive factor. Defaults to 0.1.
+
+    Returns:
+        np.ndarray: rescaled image of same size as input
+    """
+    return (mult * image) + 0.1
diff --git a/pylst/visualization/__init__.py b/pylst/visualization/__init__.py
new file mode 100644
index 0000000..ab62522
--- /dev/null
+++ b/pylst/visualization/__init__.py
@@ -0,0 +1,8 @@
+# pylst/visualization/__init__.py
+
+#from .plotter import plot_data
+from pylst.visualization.plotter import plot_data
+from pylst.visualization.plotter import show
+from pylst.visualization.plotter import rastshow
+from .plotter import create_interactive_map, histogram_equalization
+
diff --git a/pylst/visualization/plotter.py b/pylst/visualization/plotter.py
new file mode 100644
index 0000000..b834a88
--- /dev/null
+++ b/pylst/visualization/plotter.py
@@ -0,0 +1,60 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Azad Rasul
+"""
+# pylst/visualization/plotter.py
+import matplotlib.pyplot as plt
+import cv2
+import rasterio
+from rasterio.plot import show
+
+def plot_data(data):
+    # Function to plot data using matplotlib
+    plt.plot(data)
+    plt.show()
+
+from localtileserver import TileClient, get_leaflet_tile_layer
+from ipyleaflet import Map
+
+def create_interactive_map(tif_path):
+    # Create a TileClient and get the tile layer
+    client = TileClient(tif_path)
+    layer = get_leaflet_tile_layer(client)
+
+    # Create a map centered around the data
+    m = Map(center=client.center(), zoom=client.default_zoom)
+    
+    # Add the tile layer to the map
+    m.add(layer)
+    
+    return m
+
+
+
+def rastshow(path_to_raster, title='', colorbar_label='', figsize='', xlabel='', ylabel='', cmap='viridis', extent=None):
+    # Open the raster dataset
+    with rasterio.open(path_to_raster) as src:
+        # Display the raster data
+        fig, ax = plt.subplots(1, 1, figsize=figsize)
+        show(src, ax=ax, title=title, cmap=cmap, extent=extent)
+        plt.colorbar(ax=ax, mappable=ax.images[0], label=colorbar_label)
+        ax.set_xlabel(xlabel)
+        ax.set_ylabel(ylabel)
+        plt.show()
+
+def histogram_equalization(image_path):
+    # Function for histogram equalization using OpenCV
+    img = cv2.imread(image_path, 0)
+    equ = cv2.equalizeHist(img)
+
+    # Display original and equalized images side by side
+    plt.subplot(1, 2, 1)
+    plt.imshow(img, cmap='gray')
+    plt.title('Original Image')
+
+    plt.subplot(1, 2, 2)
+    plt.imshow(equ, cmap='gray')
+    plt.title('Equalized Image')
+
+    plt.show()
+