diff --git a/city_metrix/layers/impervious_surface.py b/city_metrix/layers/impervious_surface.py
index 74e624d..ea6772e 100644
--- a/city_metrix/layers/impervious_surface.py
+++ b/city_metrix/layers/impervious_surface.py
@@ -7,8 +7,14 @@
 
 
 class ImperviousSurface(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 impervious_surface
@@ -19,5 +25,5 @@ def get_data(self, bbox):
                                          .sum()
                                          )
 
-        data = get_image_collection(imperv_surf, bbox, 100, "imperv surf")
+        data = get_image_collection(imperv_surf, bbox, self.spatial_resolution, "imperv surf")
         return data.change_year_index
diff --git a/city_metrix/layers/ndvi_sentinel2_gee.py b/city_metrix/layers/ndvi_sentinel2_gee.py
index e49b3d3..abea8c1 100644
--- a/city_metrix/layers/ndvi_sentinel2_gee.py
+++ b/city_metrix/layers/ndvi_sentinel2_gee.py
@@ -7,15 +7,17 @@ class NdviSentinel2(Layer):
     Attributes:
         year: The satellite imaging year.
         spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+        ndvi_threshold: minimum NDVI value for retrieval
     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, spatial_resolution=10, **kwargs):
+    def __init__(self, year=None, spatial_resolution=10, ndvi_threshold=None, **kwargs):
         super().__init__(**kwargs)
         self.year = year
         self.spatial_resolution = spatial_resolution
+        self.ndvi_threshold = ndvi_threshold
 
     def get_data(self, bbox):
         if self.year is None:
@@ -45,4 +47,7 @@ def calculate_ndvi(image):
         ndvi_data = (get_image_collection(ic, bbox, self.spatial_resolution, "NDVI")
                      .NDVI)
 
-        return ndvi_data
+        if self.ndvi_threshold is not None:
+            return ndvi_data.where(ndvi_data >= self.ndvi_threshold)
+        else:
+            return ndvi_data
diff --git a/city_metrix/metrics/built_land_with_vegetation.py b/city_metrix/metrics/built_land_with_vegetation.py
new file mode 100644
index 0000000..3a967a9
--- /dev/null
+++ b/city_metrix/metrics/built_land_with_vegetation.py
@@ -0,0 +1,20 @@
+from geopandas import GeoDataFrame, GeoSeries
+from city_metrix.layers import EsaWorldCoverClass, EsaWorldCover, NdviSentinel2
+
+def built_land_with_vegetation(zones: GeoDataFrame) -> GeoSeries:
+    """
+    Get percentage of built up land (using ESA world cover) with NDVI vegetation cover.
+    :param zones: GeoDataFrame with geometries to collect zonal stats on
+    :return: Pandas Series of percentages
+    """
+    ndvi_sampling_year = 2020
+    ndvi_threshold = 0.4
+    built_up_land = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)
+    vegetation = NdviSentinel2(year=ndvi_sampling_year, ndvi_threshold=ndvi_threshold)
+
+    built_land = built_up_land.groupby(zones).count()
+    vegetation_cover_in_built_land = vegetation.mask(built_up_land).groupby(zones).count()
+
+    fraction_vegetation_in_built_up_land = (vegetation_cover_in_built_land.fillna(0) / built_land)
+
+    return fraction_vegetation_in_built_up_land
diff --git a/tests/conftest.py b/tests/conftest.py
index b56318c..bd3cd7d 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,34 +1,6 @@
 from city_metrix.layers import Layer
-
-import shapely.geometry as geometry
-import geopandas as gpd
 from geocube.api.core import make_geocube
-
-
-def create_fishnet_grid(min_x, min_y, max_x, max_y, cell_size):
-    x, y = (min_x, min_y)
-    geom_array = []
-
-    # Polygon Size
-    while y < (max_y - 0.000001):
-        while x < (max_x - 0.000001):
-            geom = geometry.Polygon(
-                [
-                    (x, y),
-                    (x, y + cell_size),
-                    (x + cell_size, y + cell_size),
-                    (x + cell_size, y),
-                    (x, y),
-                ]
-            )
-            geom_array.append(geom)
-            x += cell_size
-        x = min_x
-        y += cell_size
-
-    fishnet = gpd.GeoDataFrame(geom_array, columns=["geometry"]).set_crs("EPSG:4326")
-    return fishnet
-
+from tests.tools.general_tools import create_fishnet_grid
 
 # Test zones of a regular 0.01x0.01 grid over a 0.1x0.1 extent
 ZONES = create_fishnet_grid(106.7, -6.3, 106.8, -6.2, 0.01).reset_index()
diff --git a/tests/resources/bbox_constants.py b/tests/resources/bbox_constants.py
index 789ab48..331f9e4 100644
--- a/tests/resources/bbox_constants.py
+++ b/tests/resources/bbox_constants.py
@@ -2,10 +2,13 @@
 
 
 BBOX_BRA_LAURO_DE_FREITAS_1 = (
-    -38.35530428121955,
-    -12.821710300686393,
-    -38.33813814352424,
-    -12.80363249765361,
+    -38.35530428121955,-12.821710300686393,
+    -38.33813814352424,-12.80363249765361
+)
+
+BBOX_BRA_LAURO_DE_FREITAS_2 = (
+    -38.36184,-12.84134,
+    -38.34614,-12.82917
 )
 
 BBOX_BRA_SALVADOR_ADM4 = (
@@ -13,6 +16,7 @@
     -38.3041637148564007,-12.75607703449720631
 )
 
+
 # UTM Zones 22S and 23S
 BBOX_BRA_BRASILIA = (
     -48.07651,-15.89788
@@ -24,3 +28,8 @@
     -38.39993,-12.93239
 )
 
+BBOX_JAKARTA_ZONE_TEST = (
+    106, -7,
+    107, -6
+)
+
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 7ca4581..5cdd525 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
@@ -4,7 +4,7 @@
 import shutil
 from collections import namedtuple
 
-from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1
+from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1, BBOX_BRA_LAURO_DE_FREITAS_2
 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
@@ -18,23 +18,24 @@
 
 # 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
+BBOX = BBOX_BRA_LAURO_DE_FREITAS_2
 
 # 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'
+    qgis_project_file = 'layers_for_br_lauro_de_freitas2.qgz'
 
-    source_folder = os.path.dirname(__file__)
-    target_folder = get_target_folder_path()
-    create_target_folder(target_folder, True)
+    if RUN_DUMPS is True:
+        source_folder = os.path.dirname(__file__)
+        target_folder = get_target_folder_path()
+        create_target_folder(target_folder, True)
 
-    source_qgis_file = os.path.join(source_folder, qgis_project_file)
-    target_qgis_file = os.path.join(target_folder, qgis_project_file)
-    shutil.copyfile(source_qgis_file, target_qgis_file)
+        source_qgis_file = os.path.join(source_folder, qgis_project_file)
+        target_qgis_file = os.path.join(target_folder, qgis_project_file)
+        shutil.copyfile(source_qgis_file, target_qgis_file)
 
-    print("\n\033[93m QGIS project file and layer files written to folder %s.\033[0m\n" % target_folder)
+        print("\n\033[93m QGIS project file and layer files written to folder %s.\033[0m\n" % target_folder)
 
 @pytest.fixture
 def target_folder():
diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz b/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz
deleted file mode 100644
index 759515e..0000000
Binary files a/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz and /dev/null differ
diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas2.qgz b/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas2.qgz
new file mode 100644
index 0000000..57b4f5f
Binary files /dev/null and b/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas2.qgz differ
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 5e89328..e99d75c 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
@@ -1,7 +1,10 @@
 # This code is mostly intended for manual execution
 # Execution configuration is specified in the conftest file
 import pytest
+import os
 
+from city_metrix import built_land_with_high_land_surface_temperature, built_land_with_low_surface_reflectivity, \
+    built_land_without_tree_cover, natural_areas, mean_tree_cover, urban_open_space
 from city_metrix.layers import (
     Albedo,
     AlosDSM,
@@ -21,141 +24,216 @@
     TreeCanopyHeight,
     TreeCover,
     UrbanLandUse,
-    WorldPop, Layer
+    WorldPop, Layer, ImperviousSurface
 )
+from city_metrix.metrics.built_land_with_vegetation import built_land_with_vegetation
 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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'albedo.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'alos_dsm.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'average_net_building_height.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'esa_world_cover.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'high_land_surface_temperature.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'land_surface_temperature.tif')
-    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, 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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'nasa_dem.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'natural_areas.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'ndvi_sentinel2_gee.tif')
-    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, 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, 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, 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, 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, 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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'tree_canopy_height.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'tree_cover.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'urban_land_use.tif')
-    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, target_spatial_resolution_multiplier):
-    file_path = prep_output_path(target_folder, 'world_pop.tif')
-    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)
+from ...tools.general_tools import get_class_default_spatial_resolution, create_fishnet_grid, write_metric, \
+    get_zones_from_bbox_info
+
+
+class TestLayerWrite:
+    @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
+    def test_write_albedo(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'albedo.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'alos_dsm.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'average_net_building_height.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'esa_world_cover.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'high_land_surface_temperature.tif')
+        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_impervious_surface(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'impervious_surface.tif')
+        target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(ImperviousSurface())
+        LandSurfaceTemperature(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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'land_surface_temperature.tif')
+        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(self, 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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'nasa_dem.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'natural_areas.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'ndvi_sentinel2_gee.tif')
+        target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(NdviSentinel2())
+        NdviSentinel2(year=2023, ndvi_threshold =0.4, 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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'open_buildings.geojson')
+        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(self, 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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'overture_buildings.geojson')
+        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(self, 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(self, 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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'tree_canopy_height.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'tree_cover.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'urban_land_use.tif')
+        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(self, target_folder, bbox_info, target_spatial_resolution_multiplier):
+        file_path = prep_output_path(target_folder, 'world_pop.tif')
+        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)
+
+class TestMetricWrite:
+    CELL_SIZE = 0.0005
+    def test_write_built_land_with_high_land_surface_temperature(self, target_folder, bbox_info):
+        zones = get_zones_from_bbox_info(bbox_info, self.CELL_SIZE)
+        indicator = built_land_with_high_land_surface_temperature(zones)
+
+        file_path = os.path.join(target_folder,'metric_built_land_with_high_land_surface_temperature.geojson')
+        write_metric(indicator, zones, 'count', file_path)
+
+        assert verify_file_is_populated(file_path)
+
+    def test_write_built_land_with_low_surface_reflectivity(self, target_folder, bbox_info):
+        zones = get_zones_from_bbox_info(bbox_info, self.CELL_SIZE)
+        indicator = built_land_with_low_surface_reflectivity(zones)
+
+        file_path = os.path.join(target_folder,'metric_built_land_with_low_surface_reflectivity.geojson')
+        write_metric(indicator, zones, 'count', file_path)
+
+        assert verify_file_is_populated(file_path)
+
+    def test_write_built_land_with_vegetation(self, target_folder, bbox_info):
+        zones = get_zones_from_bbox_info(bbox_info, self.CELL_SIZE)
+        indicator = built_land_with_vegetation(zones)
+
+        file_path = os.path.join(target_folder,'metric_built_land_with_vegetation.geojson')
+        write_metric(indicator, zones, 'count', file_path)
+
+        assert verify_file_is_populated(file_path)
+
+    def test_write_built_land_without_tree_cover(self, target_folder, bbox_info):
+        zones = get_zones_from_bbox_info(bbox_info, self.CELL_SIZE)
+        indicator = built_land_without_tree_cover(zones)
+
+        file_path = os.path.join(target_folder,'metric_built_land_without_tree_cover.geojson')
+        write_metric(indicator, zones, 'count', file_path)
+
+        assert verify_file_is_populated(file_path)
+
+    def test_write_mean_tree_cover(self, target_folder, bbox_info):
+        zones = get_zones_from_bbox_info(bbox_info, self.CELL_SIZE)
+        indicator = mean_tree_cover(zones)
+
+        file_path = os.path.join(target_folder,'metric_mean_tree_cover.geojson')
+        write_metric(indicator, zones, 'mean', file_path)
+
+        assert verify_file_is_populated(file_path)
+
+    def test_write_natural_areas(self, target_folder, bbox_info):
+        zones = get_zones_from_bbox_info(bbox_info, self.CELL_SIZE)
+        indicator = natural_areas(zones)
+
+        file_path = os.path.join(target_folder,'metric_natural_areas.geojson')
+        write_metric(indicator, zones, 'mean', file_path)
+
+        assert verify_file_is_populated(file_path)
+
+    def test_write_urban_open_space(self, target_folder, bbox_info):
+        zones = get_zones_from_bbox_info(bbox_info, self.CELL_SIZE)
+        indicator = urban_open_space(zones)
+
+        file_path = os.path.join(target_folder,'metric_urban_open_space.geojson')
+        write_metric(indicator, zones, 'count', file_path)
+
+        assert verify_file_is_populated(file_path)
\ No newline at end of file
diff --git a/tests/test_layer_dimensions.py b/tests/test_layer_dimensions.py
index c3fdd6d..40186e3 100644
--- a/tests/test_layer_dimensions.py
+++ b/tests/test_layer_dimensions.py
@@ -1,19 +1,128 @@
-from city_metrix.layers import NdviSentinel2
+import ee
+import pytest
+
+from city_metrix.layers import NdviSentinel2, TreeCover, Albedo, AlosDSM
 from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1
+from city_metrix.layers.layer import get_image_collection
 from tests.tools.general_tools import post_process_layer
 
+EE_IMAGE_DIMENSION_TOLERANCE = 1  # Tolerance compensates for variable results from GEE service
 COUNTRY_CODE_FOR_BBOX = 'BRA'
 BBOX = BBOX_BRA_LAURO_DE_FREITAS_1
 
+def test_read_image_collection():
+    ic = ee.ImageCollection("ESA/WorldCover/v100")
+    data = get_image_collection(ic, BBOX, 10, "test")
+
+    expected_crs = 32724
+    expected_x_dimension = 187
+    expected_y_dimension = 199
+
+    assert data.rio.crs == expected_crs
+    assert (
+        pytest.approx(expected_x_dimension, rel=EE_IMAGE_DIMENSION_TOLERANCE) == "x",
+        pytest.approx(expected_y_dimension, rel=EE_IMAGE_DIMENSION_TOLERANCE) == "y"
+    )
+
+def test_read_image_collection_scale():
+    ic = ee.ImageCollection("ESA/WorldCover/v100")
+    data = get_image_collection(ic, BBOX, 100, "test")
+    expected_x_dimension = 19
+    expected_y_dimension = 20
+    assert data.dims == {"x": expected_x_dimension, "y": expected_y_dimension}
+
+def test_albedo_dimensions():
+    data = Albedo().get_data(BBOX)
+    analysis_data = post_process_layer(data, value_threshold=0.1, convert_to_percentage=True)
+
+    expected_min = 0
+    expected_max = 34
+    expected_peak_value = 15
+    # peak_value, peak_count = get_count_by_value(analysis_data, expected_min, expected_max)
+
+    # Bounding values
+    actual_min = analysis_data.values.min()
+    actual_max = analysis_data.values.max()
+
+    # Peak frequency
+    full_count = analysis_data.size
+    mid_count_pct = get_value_percent(analysis_data, expected_peak_value, full_count, 0)
+
+    # Value range
+    assert actual_min == expected_min
+    assert actual_max == expected_max
+    # Peak frequency
+    assert mid_count_pct == 21
+
+def test_alos_dsm_dimensions():
+    analysis_data = AlosDSM().get_data(BBOX)
+
+    expected_min = 16
+    expected_max = 86
+    expected_peak_value = 56
+    peak_value, peak_count = get_count_by_value(analysis_data, expected_min, expected_max)
+
+    # Bounding values
+    actual_min = analysis_data.values.min()
+    actual_max = analysis_data.values.max()
+
+    # Peak frequency
+    full_count = analysis_data.size
+    mid_count_pct = get_value_percent(analysis_data, expected_peak_value, full_count, 0)
+
+    # Value range
+    assert actual_min == expected_min
+    assert actual_max == expected_max
+    # Peak frequency
+    assert mid_count_pct == 3
+    
 def test_ndvi_dimensions():
     data = NdviSentinel2(year=2023).get_data(BBOX)
-    data_for_map = post_process_layer(data, value_threshold=0.4, convert_to_percentage=True)
+    analysis_data = post_process_layer(data, value_threshold=0.4, convert_to_percentage=True)
 
     expected_min = 0
-    actual_min = data_for_map.values.min()
     expected_max = 85
-    actual_max = data_for_map.values.max()
+    expected_peak_value = 78
+    # peak_value, peak_count = get_count_by_value(analysis_data, expected_min, expected_max)
+
+    # Bounding values
+    actual_min = analysis_data.values.min()
+    actual_max = analysis_data.values.max()
 
+    # Peak frequency
+    full_count = analysis_data.size
+    mid_count_pct = get_value_percent(analysis_data, expected_peak_value, full_count, 0)
+
+    # Value range
     assert actual_min == expected_min
     assert actual_max == expected_max
+    # Peak frequency
+    assert mid_count_pct == 11
+
+
+def test_tree_cover():
+    actual = TreeCover().get_data(BBOX).mean()
+    expected = 54.0
+    tolerance = 0.1
+    assert (
+            pytest.approx(expected, rel=tolerance) == actual
+    )
+
+def get_value_percent(data, value, full_count, precision):
+    count_for_value = data.values[data.values == value].size
+    percent_of_cells_with_value = get_rounded_pct(full_count, count_for_value, precision)
+    return percent_of_cells_with_value
+
+def get_rounded_pct(full_count, this_count, precision):
+    return round((this_count/full_count)*100, precision)
+
+def get_count_by_value(data, min_value, max_value):
+    peak_value = None
+    peak_count = 0
+    for x in range(min_value, max_value):
+        count = data.values[data.values == x].size
+        if count > peak_count:
+            peak_count = count
+            peak_value = x
 
+    return peak_value, peak_count
\ No newline at end of file
diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py
index 90942f7..ad518ee 100644
--- a/tests/test_layer_parameters.py
+++ b/tests/test_layer_parameters.py
@@ -2,8 +2,6 @@
 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,
@@ -11,198 +9,309 @@
     AverageNetBuildingHeight,
     BuiltUpHeight,
     EsaWorldCover, EsaWorldCoverClass,
+    HighLandSurfaceTemperature,
+    ImperviousSurface,
     LandSurfaceTemperature,
     NasaDEM,
     NaturalAreas,
     NdviSentinel2,
+    OpenStreetMap,
     TreeCanopyHeight,
     TreeCover,
     UrbanLandUse,
-    WorldPop, OpenStreetMap, HighLandSurfaceTemperature
+    WorldPop, OpenBuildings
 )
 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
+RESOLUTION_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
+class TestSpatialResolution:
+    """
+    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.
+    """
+
+    def test_albedo_downsampled_spatial_resolution(self):
+        class_instance = Albedo()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_alos_dsm_downsampled_spatial_resolution(self):
+        class_instance = AlosDSM()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_average_net_building_height_downsampled_spatial_resolution(self):
+        class_instance = AverageNetBuildingHeight()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_built_up_height_downsampled_spatial_resolution(self):
+        class_instance = BuiltUpHeight()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_esa_world_cover_downsampled_spatial_resolution(self):
+        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, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_high_land_surface_temperature_downsampled_spatial_resolution(self):
+        class_instance = HighLandSurfaceTemperature()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_impervious_surface_downsampled_spatial_resolution(self):
+        class_instance = ImperviousSurface()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_land_surface_temperature_downsampled_spatial_resolution(self):
+        class_instance = LandSurfaceTemperature()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_nasa_dem_downsampled_spatial_resolution(self):
+        class_instance = NasaDEM()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_natural_areas_downsampled_spatial_resolution(self):
+        class_instance = NaturalAreas()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_ndvi_sentinel2_downsampled_spatial_resolution(self):
+        class_instance = NdviSentinel2(year=2023)
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_tree_canopy_height_downsampled_spatial_resolution(self):
+        class_instance = TreeCanopyHeight()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_tree_cover_downsampled_spatial_resolution(self):
+        class_instance = TreeCover()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_urban_land_use_downsampled_spatial_resolution(self):
+        class_instance = UrbanLandUse()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_world_pop_downsampled_spatial_resolution(self):
+        class_instance = WorldPop()
+        default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = (
+            _get_sample_data(class_instance, BBOX, DOWNSIZE_FACTOR))
+        downsizing_is_within_tolerances = _evaluate_raster_value(default_res_data, downsized_res_data)
+
+        assert pytest.approx(target_downsized_res, rel=RESOLUTION_TOLERANCE) == estimated_actual_res
+        assert downsizing_is_within_tolerances
+
+    def test_halved_up_sampled_spatial_resolution(self):
+        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, BBOX)
+        estimated_actual_resolution = _estimate_spatial_resolution(modified_data)
+
+        assert expected_resolution == estimated_actual_resolution
+
+    def test_null_spatial_resolution(self):
+        class_instance = Albedo()
+        spatial_resolution=None
+
+        with pytest.raises(Exception) as e_info:
+            _get_modified_resolution_data(class_instance, spatial_resolution, BBOX)
+        msg_correct = False if str(e_info.value).find("Spatial_resolution cannot be None.") == -1 else True
+        assert msg_correct
+
+class TestOtherParameters:
+    def test_albedo_threshold(self):
+        threshold=0.1
+        data = Albedo(threshold=threshold).get_data(BBOX)
+        max_albedo = data.values[~np.isnan(data.values)].max()
+
+        assert threshold > max_albedo,\
+            f"Maximum value ({max_albedo}) in Albedo dataset is not < threshold of {threshold}."
+
+    def test_albedo_dates(self):
+        with pytest.raises(Exception) as e_info:
+            Albedo(start_date="2021-01-01", end_date="2021-01-02").get_data(BBOX)
+        msg_correct = False if str(e_info.value).find("max() arg is an empty sequence") == -1 else True
+        assert msg_correct
+
+        with pytest.raises(Exception) as e_info:
+            Albedo(start_date="2021-01-01", end_date=None).get_data(BBOX)
+        msg_correct = False if str(e_info.value).find("max() arg is an empty sequence") == -1 else True
+        assert msg_correct
+
+    def test_high_land_surface_temperature_dates(self):
+        with pytest.raises(Exception) as e_info:
+            HighLandSurfaceTemperature(start_date="2021-01-01", end_date="2021-01-02").get_data(BBOX)
+        msg_correct = False if str(e_info.value).find("Dictionary does not contain key: 'date'") == -1 else True
+        assert msg_correct
+
+    def test_land_surface_temperature_dates(self):
+        with pytest.raises(Exception) as e_info:
+            LandSurfaceTemperature(start_date="2021-01-01", end_date="2021-01-02").get_data(BBOX)
+        msg_correct = False if str(e_info.value).find("max() arg is an empty sequence") == -1 else True
+        assert msg_correct
+
+        with pytest.raises(Exception) as e_info:
+            LandSurfaceTemperature(start_date="2021-01-01", end_date=None).get_data(BBOX)
+        msg_correct = False if str(e_info.value).find("max() arg is an empty sequence") == -1 else True
+        assert msg_correct
+
+    def test_ndvi_sentinel2_dates(self):
+        with pytest.raises(Exception) as e_info:
+            data = NdviSentinel2(year=None).get_data(BBOX)
+        msg_correct = False if str(e_info.value).find("get_data() requires a year value") == -1 else True
+        assert msg_correct
+
+        with pytest.raises(Exception) as e_info:
+            NdviSentinel2(year="1970").get_data(BBOX)
+        msg_correct = False if str(e_info.value).find("max() arg is an empty sequence") == -1 else True
+        assert msg_correct
+
+    def test_open_buildings_country(self):
+        with pytest.raises(Exception) as e_info:
+            OpenBuildings(country="ZZZ").get_data(BBOX)
+
+        msg_correct = False if str(e_info.value).find("projects/sat-io/open-datasets/VIDA_COMBINED/ZZZ' not found.") == -1 else True
+        assert msg_correct
+
+    def test_tree_cover_min_max_cover(self):
+        data = TreeCover(min_tree_cover = 150).get_data(BBOX)
+        non_null_cells = data.values[~np.isnan(data)].size
+        assert non_null_cells == 0
 
-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)
+        data = TreeCover(max_tree_cover = -1).get_data(BBOX)
+        non_null_cells = data.values[~np.isnan(data)].size
+        assert non_null_cells == 0
+
+
+def test_albedo_threshold():
+    threshold=0.1
+    data = Albedo(threshold=threshold).get_data(BBOX)
+    max_albedo = data.values[~np.isnan(data.values)].max()
+
+    assert threshold > max_albedo,\
+        f"Maximum value ({max_albedo}) in Albedo dataset is not < threshold of {threshold}."
+
+def test_ndvi_sentinel2_threshold():
+    ndvi_threshold=0.4
+    data = NdviSentinel2(year=2023, ndvi_threshold=ndvi_threshold).get_data(BBOX)
+    min_ndvi = data.values[~np.isnan(data.values)].min()
 
-    assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res
-    assert downsizing_is_within_tolerances
+    assert ndvi_threshold <= min_ndvi,\
+        f"Minimum value ({min_ndvi}) in NDVI dataset is not >= threshold of {ndvi_threshold}."
 
-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')
+    is_valid, except_str = _validate_layer_instance('t')
     assert is_valid is False
-    is_valid, except_str = validate_layer_instance(EsaWorldCoverClass.BUILT_UP)
+    is_valid, except_str = _validate_layer_instance(EsaWorldCoverClass.BUILT_UP)
     assert is_valid is False
-    is_valid, except_str = validate_layer_instance(OpenStreetMap())
+    is_valid, except_str = _validate_layer_instance(OpenStreetMap())
     assert is_valid is False
-    is_valid, except_str = validate_layer_instance(Albedo(spatial_resolution = 2))
+    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)
+def _get_sample_data(class_instance, bbox, downsize_factor):
+    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_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)
+    downsized_res_data = _get_modified_resolution_data(class_instance, downsized_resolution, bbox)
+    default_res_data = _get_modified_resolution_data(class_instance, default_res, bbox)
 
-    estimated_actual_resolution = estimate_spatial_resolution(downsized_res_data)
+    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):
+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):
+def _get_modified_resolution_data(class_instance, spatial_resolution, bbox):
     class_instance.spatial_resolution = spatial_resolution
-    data = class_instance.get_data(BBOX)
+    data = class_instance.get_data(bbox)
     return data
 
-def validate_layer_instance(obj):
+def _validate_layer_instance(obj):
     is_valid = True
     except_str = None
 
@@ -224,12 +333,12 @@ def validate_layer_instance(obj):
 
     return is_valid, except_str
 
-def estimate_spatial_resolution(data):
+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 = _get_crs_from_image_data(data)
     crs_unit = crs.axis_info[0].unit_name
 
     if crs_unit == 'metre':
@@ -241,23 +350,23 @@ def estimate_spatial_resolution(data):
 
     return estimated_actual_resolution
 
-def get_populate_ratio(dataset):
+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):
+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
+
+    populated_raw_data_ratio = _get_populate_ratio(raw_data)
+    populated_downsized_data_ratio = _get_populate_ratio(raw_data)
+    ratio_diff = abs(populated_raw_data_ratio - populated_downsized_data_ratio)
+    ratio_eval = True if ratio_diff <= ratio_tolerance else False
 
     filled_raw_data = raw_data.fillna(0)
     filled_downsized_data = downsized_data.fillna(0)
@@ -284,9 +393,9 @@ def evaluate_raster_value(raw_data, downsized_data):
     matching_rmse = True if normalized_rmse < normalized_rmse_tolerance else False
 
     # Calculate and evaluate Structural Similarity Index (SSIM)
+    ssim_index_tolerance = 0.6 if (processed_downsized_data_np.size > 100 and ratio_diff <= 0.1) else 0.4
     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/test_layers.py b/tests/test_layers.py
index 27cb7c5..d26a917 100644
--- a/tests/test_layers.py
+++ b/tests/test_layers.py
@@ -1,11 +1,10 @@
-import ee
 import pytest
-
+import numpy as np
 from city_metrix.layers import (
     Albedo,
     AlosDSM,
     AverageNetBuildingHeight,
-    NdviSentinel2,
+    BuiltUpHeight,
     EsaWorldCover,
     EsaWorldCoverClass,
     HighLandSurfaceTemperature,
@@ -14,6 +13,7 @@
     LandSurfaceTemperature,
     NasaDEM,
     NaturalAreas,
+    NdviSentinel2,
     OpenBuildings,
     OpenStreetMap,
     OpenStreetMapClass,
@@ -25,143 +25,97 @@
     UrbanLandUse,
     WorldPop
 )
-from city_metrix.layers.layer import get_image_collection
 from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1
 
-EE_IMAGE_DIMENSION_TOLERANCE = 1  # Tolerance compensates for variable results from GEE service
 # Tests are implemented for the same bounding box in Brazil.
 COUNTRY_CODE_FOR_BBOX = 'BRA'
 BBOX = BBOX_BRA_LAURO_DE_FREITAS_1
 
 def test_albedo():
-    assert Albedo().get_data(BBOX).mean()
-
+    data = Albedo().get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_alos_dsm():
-    mean = AlosDSM().get_data(BBOX).mean()
-    assert mean
-
+    data = AlosDSM().get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_average_net_building_height():
-    assert AverageNetBuildingHeight().get_data(BBOX).mean()
+    data = AverageNetBuildingHeight().get_data(BBOX)
+    assert np.size(data) > 0
 
+def test_built_up_height():
+    data = BuiltUpHeight().get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_esa_world_cover():
-    count = (
-        EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)
-        .get_data(BBOX)
-        .count()
-    )
-    assert count
-
-
-def test_read_image_collection():
-    ic = ee.ImageCollection("ESA/WorldCover/v100")
-    data = get_image_collection(ic, BBOX, 10, "test")
-
-    expected_crs = 32724
-    expected_x_dimension = 187
-    expected_y_dimension = 199
-
-    assert data.rio.crs == expected_crs
-    assert (
-        pytest.approx(expected_x_dimension, rel=EE_IMAGE_DIMENSION_TOLERANCE) == "x",
-        pytest.approx(expected_y_dimension, rel=EE_IMAGE_DIMENSION_TOLERANCE) == "y"
-    )
-
-
-def test_read_image_collection_scale():
-    ic = ee.ImageCollection("ESA/WorldCover/v100")
-    data = get_image_collection(ic, BBOX, 100, "test")
-    expected_x_dimension = 19
-    expected_y_dimension = 20
-    assert data.dims == {"x": expected_x_dimension, "y": expected_y_dimension}
-
+    land_cover_class = EsaWorldCoverClass.BUILT_UP
+    data = EsaWorldCover(land_cover_class=land_cover_class).get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_high_land_surface_temperature():
     data = HighLandSurfaceTemperature().get_data(BBOX)
-    assert data.any()
-
+    assert np.size(data) > 0
 
 def test_impervious_surface():
     data = ImperviousSurface().get_data(BBOX)
-    assert data.any()
-
+    assert np.size(data) > 0
 
 def test_land_surface_temperature():
-    mean_lst = LandSurfaceTemperature().get_data(BBOX).mean()
-    assert mean_lst
-
+    data = LandSurfaceTemperature().get_data(BBOX)
+    assert np.size(data) > 0
 
 @pytest.mark.skip(reason="layer is deprecated")
 def test_landsat_collection_2():
     bands = ["blue"]
     data = LandsatCollection2(bands).get_data(BBOX)
-    assert data.any()
-
+    assert np.size(data) > 0
 
 def test_nasa_dem():
-    mean = NasaDEM().get_data(BBOX).mean()
-    assert mean
-
+    data = NasaDEM().get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_natural_areas():
     data = NaturalAreas().get_data(BBOX)
-    assert data.any()
+    assert np.size(data) > 0
 
 def test_ndvi_sentinel2():
     data = NdviSentinel2(year=2023).get_data(BBOX)
-    assert data is not None
-
+    assert np.size(data) > 0
 
 def test_openbuildings():
-    count = OpenBuildings(COUNTRY_CODE_FOR_BBOX).get_data(BBOX).count().sum()
-    assert count
-
+    data = OpenBuildings(COUNTRY_CODE_FOR_BBOX).get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_open_street_map():
-    count = (
-        OpenStreetMap(osm_class=OpenStreetMapClass.ROAD)
-        .get_data(BBOX)
-        .count()
-        .sum()
-    )
-    assert count
-
+    data = OpenStreetMap(osm_class=OpenStreetMapClass.ROAD).get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_overture_buildings():
-    count = OvertureBuildings().get_data(BBOX).count().sum()
-    assert count
-
+    data = OvertureBuildings().get_data(BBOX)
+    assert np.size(data) > 0
 
 @pytest.mark.skip(reason="layer is deprecated")
 def test_sentinel_2_level2():
     sentinel_2_bands = ["green"]
     data = Sentinel2Level2(sentinel_2_bands).get_data(BBOX)
-    assert data.any()
-
+    assert np.size(data) > 0
 
 def test_smart_surface_lulc():
-    count = SmartSurfaceLULC().get_data(BBOX).count()
-    assert count
+    data = SmartSurfaceLULC().get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_tree_canopy_height():
-    count = TreeCanopyHeight().get_data(BBOX).count()
-    assert count
+    data = TreeCanopyHeight().get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_tree_cover():
-    actual = TreeCover().get_data(BBOX).mean()
-    expected = 54.0
-    tolerance = 0.1
-    assert (
-            pytest.approx(expected, rel=tolerance) == actual
-    )
-
+    data = TreeCover().get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_urban_land_use():
-    assert UrbanLandUse().get_data(BBOX).count()
-
+    data = UrbanLandUse().get_data(BBOX)
+    assert np.size(data) > 0
 
 def test_world_pop():
     data = WorldPop().get_data(BBOX)
-    assert data.any()
+    assert np.size(data) > 0
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
index eb2c440..fe2f499 100644
--- a/tests/test_metrics.py
+++ b/tests/test_metrics.py
@@ -1,18 +1,50 @@
 from city_metrix import *
+from city_metrix.metrics.built_land_with_vegetation import built_land_with_vegetation
 from .conftest import ZONES
 
 
-def test_urban_open_space():
-    indicator = urban_open_space(ZONES)
-    assert indicator.size == 100
+def test_built_land_with_high_lst():
+    indicator = built_land_with_high_land_surface_temperature(ZONES)
+    expected_zone_size = ZONES.geometry.size
+    actual_indicator_size = indicator.size
+    assert expected_zone_size == actual_indicator_size
 
 
 def test_built_land_with_low_surface_reflectivity():
     indicator = built_land_with_low_surface_reflectivity(ZONES)
-    assert indicator.size == 100
+    expected_zone_size = ZONES.geometry.size
+    actual_indicator_size = indicator.size
+    assert expected_zone_size == actual_indicator_size
 
+def test_built_land_with_vegetation():
+    indicator = built_land_with_vegetation(ZONES)
+    expected_zone_size = ZONES.geometry.size
+    actual_indicator_size = indicator.size
+    assert expected_zone_size == actual_indicator_size
 
-def test_high_lst():
-    indicator = built_land_with_high_land_surface_temperature(ZONES)
-    assert indicator.size == 100
+def test_built_land_without_tree_cover():
+    indicator = built_land_without_tree_cover(ZONES)
+    expected_zone_size = ZONES.geometry.size
+    actual_indicator_size = indicator.size
+    assert expected_zone_size == actual_indicator_size
+
+
+def test_mean_tree_cover():
+    indicator = mean_tree_cover(ZONES)
+    expected_zone_size = ZONES.geometry.size
+    actual_indicator_size = indicator.size
+    assert expected_zone_size == actual_indicator_size
 
+
+def test_natural_areas():
+    indicator = natural_areas(ZONES)
+    expected_zone_size = ZONES.geometry.size
+    actual_indicator_size = indicator.size
+    assert expected_zone_size == actual_indicator_size
+
+
+def test_urban_open_space():
+    indicator = urban_open_space(ZONES)
+    expected_zone_size = ZONES.geometry.size
+    actual_indicator_size = indicator.size
+    assert expected_zone_size == actual_indicator_size
\ No newline at end of file
diff --git a/tests/tools/general_tools.py b/tests/tools/general_tools.py
index d3e56b8..b811ddd 100644
--- a/tests/tools/general_tools.py
+++ b/tests/tools/general_tools.py
@@ -1,6 +1,10 @@
 import os
-import tempfile
+import shapely.geometry as geometry
 import numpy as np
+import geopandas as gpd
+from geopandas import GeoDataFrame
+from city_metrix.layers.layer import get_utm_zone_epsg
+
 
 def is_valid_path(path: str):
     return os.path.exists(path)
@@ -63,4 +67,47 @@ def get_class_default_spatial_resolution(obj):
 
 def get_class_from_instance(obj):
     cls = obj.__class__()
-    return cls
\ No newline at end of file
+    return cls
+
+def write_metric(metrics, zones: GeoDataFrame, source_column_name, target_path):
+    metrics = zones.merge(metrics, left_index=True, right_index=True)
+    metrics.rename(columns={source_column_name: 'metric'}, inplace=True)
+    
+    # Project to local utm
+    zone_bounds = zones.geometry[0].envelope.bounds
+    target_epsg_crs = get_utm_zone_epsg(zone_bounds).partition('EPSG:')[2]
+    projected_metrics = metrics.to_crs(epsg=target_epsg_crs, inplace=False)
+
+    projected_metrics.to_file(target_path)
+
+def create_fishnet_grid(min_x, min_y, max_x, max_y, cell_size):
+    x, y = (min_x, min_y)
+    geom_array = []
+
+    # Polygon Size
+    while y < (max_y - 0.000001):
+        while x < (max_x - 0.000001):
+            geom = geometry.Polygon(
+                [
+                    (x, y),
+                    (x, y + cell_size),
+                    (x + cell_size, y + cell_size),
+                    (x + cell_size, y),
+                    (x, y),
+                ]
+            )
+            geom_array.append(geom)
+            x += cell_size
+        x = min_x
+        y += cell_size
+
+    fishnet = gpd.GeoDataFrame(geom_array, columns=["geometry"]).set_crs("EPSG:4326")
+    return fishnet
+
+def get_zones_from_bbox_info(bbox_info, cell_size):
+    min_x = bbox_info.bounds[0]
+    min_y = bbox_info.bounds[1]
+    max_x = bbox_info.bounds[2]
+    max_y = bbox_info.bounds[3]
+    zones = create_fishnet_grid(min_x, min_y, max_x, max_y, cell_size)
+    return zones
\ No newline at end of file