diff --git a/.github/requirements.txt b/.github/requirements.txt
index aff3f1d..f5069a9 100644
--- a/.github/requirements.txt
+++ b/.github/requirements.txt
@@ -6,7 +6,7 @@ odc-stac==0.3.8
 pystac-client==0.7.5
 pytest==7.4.3
 xarray-spatial==0.3.7
-xee==0.0.3
+xee==0.0.15
 utm==0.7.0
 osmnx==1.9.3
 dask[complete]==2023.11.0
@@ -16,5 +16,6 @@ geemap==0.32.0
 pip==23.3.1
 boto3==1.34.124
 scikit-learn==1.5.1
+scikit-image==0.24.0
 overturemaps==0.6.0
 exactextract==0.2.0.dev252
diff --git a/city_metrix/layers/albedo.py b/city_metrix/layers/albedo.py
index 7bf7b11..341786b 100644
--- a/city_metrix/layers/albedo.py
+++ b/city_metrix/layers/albedo.py
@@ -4,12 +4,20 @@
 
 from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
-
 class Albedo(Layer):
-    def __init__(self, start_date="2021-01-01", end_date="2022-01-01", threshold=None, **kwargs):
+    """
+    Attributes:
+        start_date: starting date for data retrieval
+        end_date: ending date for data retrieval
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+        threshold: threshold value for filtering the retrieval
+    """
+
+    def __init__(self, start_date="2021-01-01", end_date="2022-01-01", spatial_resolution=10, threshold=None, **kwargs):
         super().__init__(**kwargs)
         self.start_date = start_date
         self.end_date = end_date
+        self.spatial_resolution = spatial_resolution
         self.threshold = threshold
 
     def get_data(self, bbox):
@@ -30,7 +38,7 @@ def mask_and_count_clouds(s2wc, geom):
             nb_cloudy_pixels = is_cloud.reduceRegion(
                 reducer=ee.Reducer.sum().unweighted(),
                 geometry=geom,
-                scale=10,
+                scale=self.spatial_resolution,
                 maxPixels=1e9
             )
             return s2wc.updateMask(is_cloud.eq(0)).set('nb_cloudy_pixels',
@@ -88,7 +96,9 @@ def calc_s2_albedo(image):
         s2_albedo = dataset.map(calc_s2_albedo)
         albedo_mean = s2_albedo.reduce(ee.Reducer.mean())
 
-        data = get_image_collection(ee.ImageCollection(albedo_mean), bbox, 10, "albedo").albedo_mean
+        data = (get_image_collection(
+            ee.ImageCollection(albedo_mean), bbox, self.spatial_resolution, "albedo")
+                .albedo_mean)
 
         if self.threshold is not None:
             return data.where(data < self.threshold)
diff --git a/city_metrix/layers/alos_dsm.py b/city_metrix/layers/alos_dsm.py
index 95187e6..70000eb 100644
--- a/city_metrix/layers/alos_dsm.py
+++ b/city_metrix/layers/alos_dsm.py
@@ -6,8 +6,14 @@
 
 
 class AlosDSM(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=30, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         dataset = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2")
@@ -16,6 +22,6 @@ def get_data(self, bbox):
                                       .select('DSM')
                                       .mean()
                                       )
-        data = get_image_collection(alos_dsm, bbox, 30, "ALOS DSM").DSM
+        data = get_image_collection(alos_dsm, bbox, self.spatial_resolution, "ALOS DSM").DSM
 
         return data
diff --git a/city_metrix/layers/average_net_building_height.py b/city_metrix/layers/average_net_building_height.py
index a5b26f0..d0f49f2 100644
--- a/city_metrix/layers/average_net_building_height.py
+++ b/city_metrix/layers/average_net_building_height.py
@@ -5,10 +5,15 @@
 
 from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
-
 class AverageNetBuildingHeight(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=100, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         # https://ghsl.jrc.ec.europa.eu/ghs_buH2023.php
@@ -18,6 +23,8 @@ def get_data(self, bbox):
         # GLOBE - ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A")
 
         anbh = ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A")
-        data = get_image_collection(ee.ImageCollection(anbh), bbox, 100, "average net building height").b1
+        data = (get_image_collection(
+            ee.ImageCollection(anbh), bbox, self.spatial_resolution, "average net building height")
+                .b1)
         
         return data
diff --git a/city_metrix/layers/built_up_height.py b/city_metrix/layers/built_up_height.py
index 05f3159..ab080f5 100644
--- a/city_metrix/layers/built_up_height.py
+++ b/city_metrix/layers/built_up_height.py
@@ -7,8 +7,14 @@
 
 
 class BuiltUpHeight(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=100, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         # Notes for Heat project:
@@ -18,6 +24,6 @@ def get_data(self, bbox):
         # ee.ImageCollection("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_R2023A")
         
         built_height = ee.Image("JRC/GHSL/P2023A/GHS_BUILT_H/2018")
-        data = get_image_collection(ee.ImageCollection(built_height), bbox, 100, "built up height")
+        data = get_image_collection(ee.ImageCollection(built_height), bbox, self.spatial_resolution, "built up height")
         return data.built_height
 
diff --git a/city_metrix/layers/esa_world_cover.py b/city_metrix/layers/esa_world_cover.py
index ed5f823..a4b0f65 100644
--- a/city_metrix/layers/esa_world_cover.py
+++ b/city_metrix/layers/esa_world_cover.py
@@ -19,32 +19,39 @@ class EsaWorldCoverClass(Enum):
     MANGROVES = 95
     MOSS_AND_LICHEN = 100
 
-
 class EsaWorldCover(Layer):
+    """
+    Attributes:
+        land_cover_class: Enum value from EsaWorldCoverClass
+        year: year used for data retrieval
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
     STAC_CATALOG_URI = "https://services.terrascope.be/stac/"
     STAC_COLLECTION_ID = "urn:eop:VITO:ESA_WorldCover_10m_2020_AWS_V1"
     STAC_ASSET_ID = "ESA_WORLDCOVER_10M_MAP"
 
-    def __init__(self, land_cover_class=None, year=2020, **kwargs):
+    def __init__(self, land_cover_class=None, year=2020, spatial_resolution=10, **kwargs):
         super().__init__(**kwargs)
         self.land_cover_class = land_cover_class
         self.year = year
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         if self.year == 2020:
             data = get_image_collection(
-                    ee.ImageCollection("ESA/WorldCover/v100"),
-                    bbox,
-                    10,
-                    "ESA world cover"
-                ).Map
+                ee.ImageCollection("ESA/WorldCover/v100"),
+                bbox,
+                self.spatial_resolution,
+                "ESA world cover"
+            ).Map
         elif self.year == 2021:
             data = get_image_collection(
-                    ee.ImageCollection("ESA/WorldCover/v200"),
-                    bbox,
-                    10,
-                    "ESA world cover"
-                ).Map
+                ee.ImageCollection("ESA/WorldCover/v200"),
+                bbox,
+                self.spatial_resolution,
+                "ESA world cover"
+            ).Map
 
         if self.land_cover_class:
             data = data.where(data == self.land_cover_class.value)
diff --git a/city_metrix/layers/high_land_surface_temperature.py b/city_metrix/layers/high_land_surface_temperature.py
index d3943e9..5faa2ae 100644
--- a/city_metrix/layers/high_land_surface_temperature.py
+++ b/city_metrix/layers/high_land_surface_temperature.py
@@ -8,19 +8,26 @@
 
 
 class HighLandSurfaceTemperature(Layer):
+    """
+    Attributes:
+        start_date: starting date for data retrieval
+        end_date: ending date for data retrieval
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
     THRESHOLD_ADD = 3
 
-    def __init__(self, start_date="2013-01-01", end_date="2023-01-01", **kwargs):
+    def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resolution=30, **kwargs):
         super().__init__(**kwargs)
         self.start_date = start_date
         self.end_date = end_date
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         hottest_date = self.get_hottest_date(bbox)
         start_date = (hottest_date - datetime.timedelta(days=45)).strftime("%Y-%m-%d")
         end_date = (hottest_date + datetime.timedelta(days=45)).strftime("%Y-%m-%d")
 
-        lst = LandSurfaceTemperature(start_date, end_date).get_data(bbox)
+        lst = LandSurfaceTemperature(start_date, end_date, self.spatial_resolution).get_data(bbox)
         lst_mean = lst.mean(dim=['x', 'y'])
         high_lst = lst.where(lst >= (lst_mean + self.THRESHOLD_ADD))
         return high_lst
diff --git a/city_metrix/layers/land_surface_temperature.py b/city_metrix/layers/land_surface_temperature.py
index 7fff632..931cb2e 100644
--- a/city_metrix/layers/land_surface_temperature.py
+++ b/city_metrix/layers/land_surface_temperature.py
@@ -5,12 +5,19 @@
 import ee
 import xarray
 
-
 class LandSurfaceTemperature(Layer):
-    def __init__(self, start_date="2013-01-01", end_date="2023-01-01", **kwargs):
+    """
+    Attributes:
+        start_date: starting date for data retrieval
+        end_date: ending date for data retrieval
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resolution=30, **kwargs):
         super().__init__(**kwargs)
         self.start_date = start_date
         self.end_date = end_date
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         def cloud_mask(image):
@@ -31,5 +38,5 @@ def apply_scale_factors(image):
             .map(apply_scale_factors) \
             .reduce(ee.Reducer.mean())
 
-        data = get_image_collection(ee.ImageCollection(l8_st), bbox, 30, "LST").ST_B10_mean
+        data = get_image_collection(ee.ImageCollection(l8_st), bbox, self.spatial_resolution, "LST").ST_B10_mean
         return data
diff --git a/city_metrix/layers/landsat_collection_2.py b/city_metrix/layers/landsat_collection_2.py
index 248227a..d82180d 100644
--- a/city_metrix/layers/landsat_collection_2.py
+++ b/city_metrix/layers/landsat_collection_2.py
@@ -29,7 +29,8 @@ def get_data(self, bbox):
             fail_on_error=False,
         )
 
-        # TODO: Determine how to output xarray
-
         qa_lst = lc2.where((lc2.qa_pixel & 24) == 0)
         return qa_lst.drop_vars("qa_pixel")
+
+
+
diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py
index 01ad6e4..c86e9fe 100644
--- a/city_metrix/layers/layer.py
+++ b/city_metrix/layers/layer.py
@@ -274,9 +274,12 @@ def get_image_collection(
     :param name: optional name to print while reporting progress
     :return:
     """
+    if scale is None:
+        raise Exception("Spatial_resolution cannot be None.")
 
     crs = get_utm_zone_epsg(bbox)
 
+    # See link regarding bug in crs specification https://github.com/google/Xee/issues/118
     ds = xr.open_dataset(
         image_collection,
         engine='ee',
diff --git a/city_metrix/layers/nasa_dem.py b/city_metrix/layers/nasa_dem.py
index 49b2d13..b5ac45d 100644
--- a/city_metrix/layers/nasa_dem.py
+++ b/city_metrix/layers/nasa_dem.py
@@ -6,8 +6,14 @@
 
 
 class NasaDEM(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=30, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         dataset = ee.Image("NASA/NASADEM_HGT/001")
@@ -16,6 +22,6 @@ def get_data(self, bbox):
                                       .select('elevation')
                                       .mean()
                                       )
-        data = get_image_collection(nasa_dem, bbox, 30, "NASA DEM").elevation
+        data = get_image_collection(nasa_dem, bbox, self.spatial_resolution, "NASA DEM").elevation
         
         return data
diff --git a/city_metrix/layers/natural_areas.py b/city_metrix/layers/natural_areas.py
index 1dfbc73..5efe4a8 100644
--- a/city_metrix/layers/natural_areas.py
+++ b/city_metrix/layers/natural_areas.py
@@ -4,13 +4,18 @@
 from .layer import Layer
 from .esa_world_cover import EsaWorldCover, EsaWorldCoverClass
 
-
 class NaturalAreas(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=10, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
-        esa_world_cover = EsaWorldCover().get_data(bbox)
+        esa_world_cover = EsaWorldCover(spatial_resolution=self.spatial_resolution).get_data(bbox)
         reclass_map = {
             EsaWorldCoverClass.TREE_COVER.value: 1,
             EsaWorldCoverClass.SHRUBLAND.value: 1,
diff --git a/city_metrix/layers/ndvi_sentinel2_gee.py b/city_metrix/layers/ndvi_sentinel2_gee.py
index c5b21b9..e49b3d3 100644
--- a/city_metrix/layers/ndvi_sentinel2_gee.py
+++ b/city_metrix/layers/ndvi_sentinel2_gee.py
@@ -4,15 +4,18 @@
 class NdviSentinel2(Layer):
     """"
     NDVI = Sentinel-2 Normalized Difference Vegetation Index
-    param: year: The satellite imaging year.
+    Attributes:
+        year: The satellite imaging year.
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
     return: a rioxarray-format DataArray
     Author of associated Jupyter notebook: Ted.Wong@wri.org
     Notebook: https://github.com/wri/cities-cities4forests-indicators/blob/dev-eric/scripts/extract-VegetationCover.ipynb
     Reference: https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index
     """
-    def __init__(self, year=None, **kwargs):
+    def __init__(self, year=None, spatial_resolution=10, **kwargs):
         super().__init__(**kwargs)
         self.year = year
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         if self.year is None:
@@ -39,8 +42,7 @@ def calculate_ndvi(image):
         ndvi_mosaic = ndvi.qualityMosaic('NDVI')
 
         ic = ee.ImageCollection(ndvi_mosaic)
-        ndvi_data = get_image_collection(ic, bbox, 10, "NDVI")
+        ndvi_data = (get_image_collection(ic, bbox, self.spatial_resolution, "NDVI")
+                     .NDVI)
 
-        xdata = ndvi_data.to_dataarray()
-
-        return xdata
+        return ndvi_data
diff --git a/city_metrix/layers/sentinel_2_level_2.py b/city_metrix/layers/sentinel_2_level_2.py
index a609293..a7ae944 100644
--- a/city_metrix/layers/sentinel_2_level_2.py
+++ b/city_metrix/layers/sentinel_2_level_2.py
@@ -50,6 +50,4 @@ def get_data(self, bbox):
         cloud_masked = s2.where(s2 != 0).where(s2.scl != 3).where(s2.scl != 8).where(s2.scl != 9).where(
             s2.scl != 10)
 
-        # TODO: Determine how to output as an xarray
-
         return cloud_masked.drop_vars("scl")
diff --git a/city_metrix/layers/tree_canopy_height.py b/city_metrix/layers/tree_canopy_height.py
index 2785272..fee1697 100644
--- a/city_metrix/layers/tree_canopy_height.py
+++ b/city_metrix/layers/tree_canopy_height.py
@@ -5,19 +5,25 @@
 import xee
 import ee
 
-
 class TreeCanopyHeight(Layer):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
     name = "tree_canopy_height"
     NO_DATA_VALUE = 0
 
-    def __init__(self, **kwargs):
+    def __init__(self, spatial_resolution=1, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         canopy_ht = ee.ImageCollection("projects/meta-forest-monitoring-okw37/assets/CanopyHeight")
         # aggregate time series into a single image
         canopy_ht = canopy_ht.reduce(ee.Reducer.mean()).rename("cover_code")
 
-        data = get_image_collection(ee.ImageCollection(canopy_ht), bbox, 1, "tree canopy height")
+        data = get_image_collection(ee.ImageCollection(canopy_ht), bbox,
+                                    self.spatial_resolution, "tree canopy height")
 
         return data.cover_code
diff --git a/city_metrix/layers/tree_cover.py b/city_metrix/layers/tree_cover.py
index fadc9bb..98bc481 100644
--- a/city_metrix/layers/tree_cover.py
+++ b/city_metrix/layers/tree_cover.py
@@ -5,18 +5,22 @@
 import xee
 import ee
 
-
 class TreeCover(Layer):
     """
     Merged tropical and nontropical tree cover from WRI
+    Attributes:
+        min_tree_cover: minimum tree-cover values used for filtering results
+        max_tree_cover: maximum tree-cover values used for filtering results
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
     """
-    
+
     NO_DATA_VALUE = 255
 
-    def __init__(self, min_tree_cover=None, max_tree_cover=None, **kwargs):
+    def __init__(self, min_tree_cover=None, max_tree_cover=None, spatial_resolution=10, **kwargs):
         super().__init__(**kwargs)
         self.min_tree_cover = min_tree_cover
         self.max_tree_cover = max_tree_cover
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         tropics = ee.ImageCollection('projects/wri-datalab/TropicalTreeCover')
@@ -24,7 +28,7 @@ def get_data(self, bbox):
         merged_ttc = tropics.merge(nontropics)
         ttc_image = merged_ttc.reduce(ee.Reducer.mean()).rename('ttc')
 
-        data = get_image_collection(ee.ImageCollection(ttc_image), bbox, 10, "tree cover").ttc
+        data = get_image_collection(ee.ImageCollection(ttc_image), bbox, self.spatial_resolution, "tree cover").ttc
 
         if self.min_tree_cover is not None:
             data = data.where(data >= self.min_tree_cover)
diff --git a/city_metrix/layers/urban_land_use.py b/city_metrix/layers/urban_land_use.py
index f34408a..fe69c75 100644
--- a/city_metrix/layers/urban_land_use.py
+++ b/city_metrix/layers/urban_land_use.py
@@ -7,9 +7,16 @@
 
 
 class UrbanLandUse(Layer):
-    def __init__(self, band='lulc', **kwargs):
+    """
+    Attributes:
+        band: raster band used for data retrieval
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, band='lulc', spatial_resolution=5, **kwargs):
         super().__init__(**kwargs)
         self.band = band
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         dataset = ee.ImageCollection("projects/wri-datalab/cities/urban_land_use/V1")
@@ -27,6 +34,6 @@ def get_data(self, bbox):
                                      .rename('lulc')
                                      )
 
-        data = get_image_collection(ulu, bbox, 5, "urban land use").lulc
+        data = get_image_collection(ulu, bbox, self.spatial_resolution, "urban land use").lulc
 
         return data
diff --git a/city_metrix/layers/world_pop.py b/city_metrix/layers/world_pop.py
index f3cf1a2..dff494a 100644
--- a/city_metrix/layers/world_pop.py
+++ b/city_metrix/layers/world_pop.py
@@ -5,10 +5,15 @@
 
 from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
-
 class WorldPop(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=100, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         # load population
@@ -20,5 +25,5 @@ def get_data(self, bbox):
                      .mean()
                      )
 
-        data = get_image_collection(world_pop, bbox, 100, "world pop")
+        data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop")
         return data.population
diff --git a/environment.yml b/environment.yml
index 37b823c..a064834 100644
--- a/environment.yml
+++ b/environment.yml
@@ -5,13 +5,13 @@ dependencies:
   - python=3.10
   - earthengine-api=0.1.379
   - geocube=0.4.2
-  - geopandas=0.14.1
+  - geopandas=0.14.4
   - rioxarray=0.15.0
   - odc-stac=0.3.8
   - pystac-client=0.7.5
   - pytest=7.4.3
   - xarray-spatial=0.3.7
-  - xee=0.0.3
+  - xee=0.0.15
   - utm=0.7.0
   - osmnx=1.9.0
   - dask[complete]=2023.11.0
@@ -22,6 +22,7 @@ dependencies:
   - pip=23.3.1
   - boto3=1.34.124
   - scikit-learn=1.5.1
+  - scikit-image==0.24.0
   - exactextract=0.2.0.dev252
   - pip:
       - overturemaps==0.6.0
diff --git a/setup.py b/setup.py
index 960c1f4..7d9242c 100644
--- a/setup.py
+++ b/setup.py
@@ -31,5 +31,6 @@
         "exactextract<=0.2.0.dev252",
         "overturemaps",
         "scikit-learn>=1.5.1",
+        "scikit-image>=0.24.0"
     ],
 )
diff --git a/tests/resources/bbox_constants.py b/tests/resources/bbox_constants.py
index 9aaf8ad..789ab48 100644
--- a/tests/resources/bbox_constants.py
+++ b/tests/resources/bbox_constants.py
@@ -9,10 +9,14 @@
 )
 
 BBOX_BRA_SALVADOR_ADM4 = (
-    -38.647320153390055,
-    -13.01748678217598787,
-    -38.3041637148564007,
-    -12.75607703449720631
+    -38.647320153390055,-13.01748678217598787,
+    -38.3041637148564007,-12.75607703449720631
+)
+
+# UTM Zones 22S and 23S
+BBOX_BRA_BRASILIA = (
+    -48.07651,-15.89788
+    -47.83736,-15.71919
 )
 
 BBOX_SMALL_TEST = (
diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py
index 5882053..7ca4581 100644
--- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py
+++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py
@@ -5,20 +5,24 @@
 from collections import namedtuple
 
 from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1
-from tools.general_tools import create_target_folder, is_valid_path
+from tests.tools.general_tools import create_target_folder, is_valid_path
 
 # RUN_DUMPS is the master control for whether the writes and tests are executed
 # Setting RUN_DUMPS to True turns on code execution.
 # Values should normally be set to False in order to avoid unnecessary execution.
 RUN_DUMPS = False
 
-# Specify None to write to a temporary default folder otherwise specify a valid custom target path.
-CUSTOM_DUMP_DIRECTORY = None
+# Multiplier applied to the default spatial_resolution of the layer
+# Use value of 1 for default resolution.
+SPATIAL_RESOLUTION_MULTIPLIER = 1
 
 # Both the tests and QGIS file are implemented for the same bounding box in Brazil.
 COUNTRY_CODE_FOR_BBOX = 'BRA'
 BBOX = BBOX_BRA_LAURO_DE_FREITAS_1
 
+# Specify None to write to a temporary default folder otherwise specify a valid custom target path.
+CUSTOM_DUMP_DIRECTORY = None
+
 def pytest_configure(config):
     qgis_project_file = 'layers_for_br_lauro_de_freitas.qgz'
 
@@ -36,6 +40,10 @@ def pytest_configure(config):
 def target_folder():
     return get_target_folder_path()
 
+@pytest.fixture
+def target_spatial_resolution_multiplier():
+    return SPATIAL_RESOLUTION_MULTIPLIER
+
 @pytest.fixture
 def bbox_info():
     bbox = namedtuple('bbox', ['bounds', 'country'])
diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py
index 5e0efb9..5e89328 100644
--- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py
+++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py
@@ -24,122 +24,138 @@
     WorldPop, Layer
 )
 from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated
+from ...tools.general_tools import get_class_default_spatial_resolution
+
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_albedo(target_folder, bbox_info):
+def test_write_albedo(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'albedo.tif')
-    Albedo().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(Albedo())
+    Albedo(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_alos_dsm(target_folder, bbox_info):
+def test_write_alos_dsm(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'alos_dsm.tif')
-    AlosDSM().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(AlosDSM())
+    AlosDSM(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_average_net_building_height(target_folder, bbox_info):
+def test_write_average_net_building_height(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'average_net_building_height.tif')
-    AverageNetBuildingHeight().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(AverageNetBuildingHeight())
+    AverageNetBuildingHeight(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_esa_world_cover(target_folder, bbox_info):
+def test_write_esa_world_cover(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'esa_world_cover.tif')
-    EsaWorldCover().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(EsaWorldCover())
+    EsaWorldCover(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_high_land_surface_temperature(target_folder, bbox_info):
+def test_write_high_land_surface_temperature(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'high_land_surface_temperature.tif')
-    HighLandSurfaceTemperature().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(HighLandSurfaceTemperature())
+    HighLandSurfaceTemperature(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_land_surface_temperature(target_folder, bbox_info):
+def test_write_land_surface_temperature(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'land_surface_temperature.tif')
-    LandSurfaceTemperature().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(LandSurfaceTemperature())
+    LandSurfaceTemperature(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 # TODO Class is no longer used, but may be useful later
 # @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-# def test_write_landsat_collection_2(target_folder, bbox_info):
+# def test_write_landsat_collection_2(target_folder, bbox_info, target_spatial_resolution_multiplier):
 #     file_path = prep_output_path(target_folder, 'landsat_collection2.tif')
 #     bands = ['green']
 #     LandsatCollection2(bands).write(bbox_info.bounds, file_path, tile_degrees=None)
 #     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_nasa_dem(target_folder, bbox_info):
+def test_write_nasa_dem(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'nasa_dem.tif')
-    NasaDEM().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(NasaDEM())
+    NasaDEM(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_natural_areas(target_folder, bbox_info):
+def test_write_natural_areas(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'natural_areas.tif')
-    NaturalAreas().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(NaturalAreas())
+    NaturalAreas(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_ndvi_sentinel2_gee(target_folder, bbox_info):
+def test_write_ndvi_sentinel2_gee(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'ndvi_sentinel2_gee.tif')
-    NdviSentinel2(year=2023).write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(NdviSentinel2())
+    NdviSentinel2(year=2023, spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_openbuildings(target_folder, bbox_info):
+def test_write_openbuildings(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'open_buildings.tif')
     OpenBuildings(bbox_info.country).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 # TODO Class write is not functional. Is class still needed or have we switched to overture?
 # @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-# def test_write_open_street_map(target_folder, bbox_info):
+# def test_write_open_street_map(target_folder, bbox_info, target_spatial_resolution_multiplier):
 #     file_path = prep_output_path(target_folder, 'open_street_map.tif')
 #     OpenStreetMap().write(bbox_info.bounds, file_path, tile_degrees=None)
 #     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_overture_buildings(target_folder, bbox_info):
+def test_write_overture_buildings(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'overture_buildings.tif')
     OvertureBuildings().write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 # TODO Class is no longer used, but may be useful later
 # @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-# def test_write_sentinel_2_level2(target_folder, bbox_info):
+# def test_write_sentinel_2_level2(target_folder, bbox_info, target_spatial_resolution_multiplier):
 #     file_path = prep_output_path(target_folder, 'sentinel_2_level2.tif')
 #     sentinel_2_bands = ["green"]
 #     Sentinel2Level2(sentinel_2_bands).write(bbox_info.bounds, file_path, tile_degrees=None)
 #     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_smart_surface_lulc(target_folder, bbox_info):
+def test_write_smart_surface_lulc(target_folder, bbox_info, target_spatial_resolution_multiplier):
+    # Note: spatial_resolution not implemented for this raster class
     file_path = prep_output_path(target_folder, 'smart_surface_lulc.tif')
     SmartSurfaceLULC().write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_tree_canopy_height(target_folder, bbox_info):
+def test_write_tree_canopy_height(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'tree_canopy_height.tif')
-    TreeCanopyHeight().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(TreeCanopyHeight())
+    TreeCanopyHeight(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_tree_cover(target_folder, bbox_info):
+def test_write_tree_cover(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'tree_cover.tif')
-    TreeCover().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(TreeCover())
+    TreeCover(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_urban_land_use(target_folder, bbox_info):
+def test_write_urban_land_use(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'urban_land_use.tif')
-    UrbanLandUse().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(UrbanLandUse())
+    UrbanLandUse(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
-def test_write_world_pop(target_folder, bbox_info):
+def test_write_world_pop(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, 'world_pop.tif')
-    WorldPop().write(bbox_info.bounds, file_path, tile_degrees=None)
+    target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(WorldPop())
+    WorldPop(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None)
     assert verify_file_is_populated(file_path)
diff --git a/tests/test_layer_dimensions.py b/tests/test_layer_dimensions.py
index 15768d6..c3fdd6d 100644
--- a/tests/test_layer_dimensions.py
+++ b/tests/test_layer_dimensions.py
@@ -1,6 +1,6 @@
 from city_metrix.layers import NdviSentinel2
 from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1
-from tests.tools import post_process_layer
+from tests.tools.general_tools import post_process_layer
 
 COUNTRY_CODE_FOR_BBOX = 'BRA'
 BBOX = BBOX_BRA_LAURO_DE_FREITAS_1
@@ -16,3 +16,4 @@ def test_ndvi_dimensions():
 
     assert actual_min == expected_min
     assert actual_max == expected_max
+
diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py
new file mode 100644
index 0000000..90942f7
--- /dev/null
+++ b/tests/test_layer_parameters.py
@@ -0,0 +1,292 @@
+import pytest
+import numpy as np
+from skimage.metrics import structural_similarity as ssim
+from pyproj import CRS
+from xrspatial import quantile
+
+from city_metrix.layers import (
+    Layer,
+    Albedo,
+    AlosDSM,
+    AverageNetBuildingHeight,
+    BuiltUpHeight,
+    EsaWorldCover, EsaWorldCoverClass,
+    LandSurfaceTemperature,
+    NasaDEM,
+    NaturalAreas,
+    NdviSentinel2,
+    TreeCanopyHeight,
+    TreeCover,
+    UrbanLandUse,
+    WorldPop, OpenStreetMap, HighLandSurfaceTemperature
+)
+from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1
+from tests.tools.general_tools import get_class_from_instance, get_class_default_spatial_resolution
+
+"""
+Evaluation of spatial_resolution property 
+To add a test for a scalable layer that has the spatial_resolution property:
+1. Add the class name to the city_metrix.layers import statement above
+2. Copy an existing test_*_spatial_resolution() test
+   a. rename for the new layer
+   b. specify a minimal class instance for the layer, not specifying the spatial_resolution attribute.
+"""
+
+COUNTRY_CODE_FOR_BBOX = 'BRA'
+BBOX = BBOX_BRA_LAURO_DE_FREITAS_1
+RESOLUTION_COMPARISON_TOLERANCE = 1
+DOWNSIZE_FACTOR = 2
+
+def test_albedo_downsampled_spatial_resolution():
+    class_instance = Albedo()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_alos_dsm_downsampled_spatial_resolution():
+    class_instance = AlosDSM()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_average_net_building_height_downsampled_spatial_resolution():
+    class_instance = AverageNetBuildingHeight()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_built_up_height_downsampled_spatial_resolution():
+    class_instance = BuiltUpHeight()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_esa_world_cover_downsampled_spatial_resolution():
+    class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_high_land_surface_temperature_downsampled_spatial_resolution():
+    class_instance = HighLandSurfaceTemperature()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_land_surface_temperature_downsampled_spatial_resolution():
+    class_instance = LandSurfaceTemperature()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_nasa_dem_downsampled_spatial_resolution():
+    class_instance = NasaDEM()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_natural_areas_downsampled_spatial_resolution():
+    class_instance = NaturalAreas()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_ndvi_sentinel2_downsampled_spatial_resolution():
+    class_instance = NdviSentinel2(year=2023)
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_tree_canopy_height_downsampled_spatial_resolution():
+    class_instance = TreeCanopyHeight()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_tree_cover_downsampled_spatial_resolution():
+    class_instance = TreeCover()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_urban_land_use_downsampled_spatial_resolution():
+    class_instance = UrbanLandUse()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_world_pop_downsampled_spatial_resolution():
+    class_instance = WorldPop()
+    default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance)
+    downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data)
+
+    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
+    assert downsizing_is_within_tolerances
+
+def test_halved_up_sampled_spatial_resolution():
+    class_instance = Albedo()
+    halved_default_resolution = get_class_default_spatial_resolution(class_instance) / 2
+
+    expected_resolution = halved_default_resolution
+    modified_data = get_modified_resolution_data(class_instance, halved_default_resolution)
+    estimated_actual_resolution = estimate_spatial_resolution(modified_data)
+
+    assert expected_resolution == estimated_actual_resolution
+
+def test_null_spatial_resolution():
+    class_instance = Albedo()
+    spatial_resolution=None
+
+    with pytest.raises(Exception) as e_info:
+        get_modified_resolution_data(class_instance, spatial_resolution)
+
+def test_function_validate_layer_instance():
+    is_valid, except_str = validate_layer_instance('t')
+    assert is_valid is False
+    is_valid, except_str = validate_layer_instance(EsaWorldCoverClass.BUILT_UP)
+    assert is_valid is False
+    is_valid, except_str = validate_layer_instance(OpenStreetMap())
+    assert is_valid is False
+    is_valid, except_str = validate_layer_instance(Albedo(spatial_resolution = 2))
+    assert is_valid is False
+
+def get_sample_data(class_instance):
+    is_valid, except_str = validate_layer_instance(class_instance)
+    if is_valid is False:
+        raise Exception(except_str)
+
+    default_res = get_class_default_spatial_resolution(class_instance)
+    downsized_resolution = DOWNSIZE_FACTOR * default_res
+
+    downsized_res_data = get_modified_resolution_data(class_instance, downsized_resolution)
+    default_res_data = get_modified_resolution_data(class_instance, default_res)
+
+    estimated_actual_resolution = estimate_spatial_resolution(downsized_res_data)
+
+    return default_res_data, downsized_res_data, downsized_resolution, estimated_actual_resolution
+
+def get_crs_from_image_data(image_data):
+    crs_string = image_data.rio.crs.data['init']
+    crs = CRS.from_string(crs_string)
+    return crs
+
+
+def get_modified_resolution_data(class_instance, spatial_resolution):
+    class_instance.spatial_resolution = spatial_resolution
+    data = class_instance.get_data(BBOX)
+    return data
+
+def validate_layer_instance(obj):
+    is_valid = True
+    except_str = None
+
+    if not obj.__class__.__bases__[0] == Layer:
+        is_valid = False
+        except_str = "Specified object '%s' is not a valid Layer class instance." % obj
+    else:
+        cls = get_class_from_instance(obj)
+        cls_name = type(cls).__name__
+        if not hasattr(obj, 'spatial_resolution'):
+            is_valid = False
+            except_str = "Class '%s' does not have spatial_resolution property." % cls_name
+        elif not obj.spatial_resolution == cls.spatial_resolution:
+            is_valid = False
+            except_str = "Do not specify spatial_resolution property value for class '%s'." % cls_name
+        elif cls.spatial_resolution is None:
+            is_valid = False
+            except_str = "Signature of class %s must specify a non-null default value for spatial_resolution. Please correct." % cls_name
+
+    return is_valid, except_str
+
+def estimate_spatial_resolution(data):
+    y_cells = float(data['y'].size - 1)
+    y_min = data.coords['y'].values.min()
+    y_max = data.coords['y'].values.max()
+
+    crs = get_crs_from_image_data(data)
+    crs_unit = crs.axis_info[0].unit_name
+
+    if crs_unit == 'metre':
+        diff_distance = y_max - y_min
+    else:
+        raise Exception('Unhandled projection units: %s for projection: %s' % (crs_unit, crs.srs))
+
+    estimated_actual_resolution = round(diff_distance / y_cells)
+
+    return estimated_actual_resolution
+
+def get_populate_ratio(dataset):
+    raw_data_size = dataset.values.size
+    populated_raw_data = dataset.values[(~np.isnan(dataset.values)) & (dataset.values > 0)]
+    populated_data_raw_size = populated_raw_data.size
+    populated_raw_data_ratio = populated_data_raw_size/raw_data_size
+    return populated_raw_data_ratio
+
+def evaluate_raster_value(raw_data, downsized_data):
+    # Below values where determined through trial and error evaluation of results in QGIS
+    ratio_tolerance = 0.2
+    normalized_rmse_tolerance = 0.3
+    ssim_index_tolerance = 0.6
+
+    populated_raw_data_ratio = get_populate_ratio(raw_data)
+    populated_downsized_data_ratio = get_populate_ratio(raw_data)
+    diff = abs(populated_raw_data_ratio - populated_downsized_data_ratio)
+    ratio_eval = True if diff <= ratio_tolerance else False
+
+    filled_raw_data = raw_data.fillna(0)
+    filled_downsized_data = downsized_data.fillna(0)
+
+    # Resample raw_data to match the resolution of the downsized_data.
+    # This operation is necessary in order to use RMSE and SSIM since they require matching array dimensions.
+    # Interpolation is set to quadratic method to intentionally not match method used by xee. (TODO Find documentation)
+    # For future releases of this method, we may need to control the interpolation to match what was used
+    # for a specific layer.
+    # Note: Initial investigations using the unresampled-raw and downsized data with evaluation by
+    # mean value, quantiles,and standard deviation were unsuccessful due to false failures on valid downsized images.
+    resampled_filled_raw_data = (filled_raw_data
+                                 .interp_like(filled_downsized_data, method='quadratic')
+                                 .fillna(0))
+
+    # Convert xarray DataArrays to numpy arrays
+    processed_raw_data_np = resampled_filled_raw_data.values
+    processed_downsized_data_np = filled_downsized_data.values
+
+    # Calculate and evaluate normalized Root Mean Squared Error (RMSE)
+    max_val = processed_downsized_data_np.max() \
+        if processed_downsized_data_np.max() > processed_raw_data_np.max() else processed_raw_data_np.max()
+    normalized_rmse = np.sqrt(np.mean(np.square(processed_downsized_data_np - processed_raw_data_np))) / max_val
+    matching_rmse = True if normalized_rmse < normalized_rmse_tolerance else False
+
+    # Calculate and evaluate Structural Similarity Index (SSIM)
+    ssim_index, _ = ssim(processed_downsized_data_np, processed_raw_data_np, full=True, data_range=max_val)
+    matching_ssim = True if ssim_index > ssim_index_tolerance else False
+
+    results_match = True if (ratio_eval & matching_rmse & matching_ssim) else False
+
+    return results_match
diff --git a/tests/tools/__init__.py b/tests/tools/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/tools.py b/tests/tools/general_tools.py
similarity index 55%
rename from tests/tools.py
rename to tests/tools/general_tools.py
index 99425df..d3e56b8 100644
--- a/tests/tools.py
+++ b/tests/tools/general_tools.py
@@ -1,5 +1,27 @@
+import os
+import tempfile
 import numpy as np
 
+def is_valid_path(path: str):
+    return os.path.exists(path)
+
+def create_target_folder(folder_path, delete_existing_files: bool):
+    if os.path.isdir(folder_path) is False:
+        os.makedirs(folder_path)
+    elif delete_existing_files is True:
+        remove_all_files_in_directory(folder_path)
+
+def remove_all_files_in_directory(directory):
+    # Iterate over all the files in the directory
+    for filename in os.listdir(directory):
+        file_path = os.path.join(directory, filename)
+        try:
+            # Check if it is a file and remove it
+            if os.path.isfile(file_path):
+                os.remove(file_path)
+        except Exception as e:
+            print(f"Error: {e}")
+
 def post_process_layer(data, value_threshold=0.4, convert_to_percentage=True):
     """
     Applies the standard post-processing adjustment used for rendering of NDVI including masking
@@ -33,3 +55,12 @@ def convert_ratio_to_percentage(data):
     values_as_percent.rio.write_crs(source_crs, inplace=True)
 
     return values_as_percent
+
+def get_class_default_spatial_resolution(obj):
+    cls = get_class_from_instance(obj)
+    default_spatial_resolution = cls.spatial_resolution
+    return default_spatial_resolution
+
+def get_class_from_instance(obj):
+    cls = obj.__class__()
+    return cls
\ No newline at end of file
diff --git a/tools/general_tools.py b/tools/general_tools.py
deleted file mode 100644
index b38d1b5..0000000
--- a/tools/general_tools.py
+++ /dev/null
@@ -1,23 +0,0 @@
-import os
-import tempfile
-
-
-def is_valid_path(path: str):
-    return os.path.exists(path)
-
-def create_target_folder(folder_path, delete_existing_files: bool):
-    if os.path.isdir(folder_path) is False:
-        os.makedirs(folder_path)
-    elif delete_existing_files is True:
-        remove_all_files_in_directory(folder_path)
-
-def remove_all_files_in_directory(directory):
-    # Iterate over all the files in the directory
-    for filename in os.listdir(directory):
-        file_path = os.path.join(directory, filename)
-        try:
-            # Check if it is a file and remove it
-            if os.path.isfile(file_path):
-                os.remove(file_path)
-        except Exception as e:
-            print(f"Error: {e}")
\ No newline at end of file