Repository: leftfield-geospatial/geedim (language: python)

Issue:
Sentinel-2 masking fails when missing cloud data.
I am having trouble using geedim with a specific image from the harmonized sentinel 2 collection, namely `COPERNICUS/S2_SR_HARMONIZED/20230809T201859_20230809T201853_T10WEE`. Whenever this images comes up in a search collection, i cannot process the collection for example to extract data with *error: Image.rename: Parameter 'input' is required.* 

I do not know whats specifically different with this image. I drilled that down to a problem with the processing geedim applies to the image when using the region parameter but i cannot find the problem itself.

i am using a conda environment with 
```
python                    3.11.9               he1021f5_0
earthengine-api           0.1.408            pyhd8ed1ab_0    conda-forge
geedim                    1.8.0              pyhd8ed1ab_0    conda-forge
```

I tested it with a bunch of images in a specific region and different region polygons. I made an example using getInfo() on the resulting image to trigger the processing, but it applies to any other client-side GEE method including geedims download():

```python
import ee
import geedim

ee.Initialize()

poly_coords1 = [
    [-121.30224921178123,71.77255565243381],
    [-121.69226386021873,71.58260352221238],
    [-120.87378241490623,71.52698680217173],
    [-120.86828925084373,71.70370209220792],
    [-121.30224921178123,71.77255565243381],
]
poly_coords2 = [
    [-121.65930487584373,71.92795376863671],
    [-121.70874335240623,71.82200371558541],
    [-121.16492011021873,71.84769170326668],
    [-121.45056464146873,71.90237555398673],
    [-121.65930487584373,71.92795376863671],
]
poly_coords3 = [
    [-120.97265936803123,71.98070472145938],
    [-120.98913886021873,71.73784463314388],
    [-120.08276678990623,71.74472884706717],
    [-120.28052069615623,71.98750044350669],
    [-120.97265936803123,71.98070472145938],
]

bounds = [
    ee.Geometry.Polygon(list(poly_coords1), proj=None, evenOdd=False),
    ee.Geometry.Polygon(list(poly_coords2), proj=None, evenOdd=False),
    ee.Geometry.Polygon(list(poly_coords3), proj=None, evenOdd=False),
]

images = [
    "20230801T200901_20230801T200857_T10WEE",
    "20230802T202849_20230802T202847_T10WEE",
    "20230803T195859_20230803T200054_T10WEE",
    "20230804T201851_20230804T201853_T10WEE",
    "20230804T201851_20230804T202447_T10WEE",
    "20230806T200859_20230806T201056_T10WEE",
    "20230807T202851_20230807T203151_T10WEE",
    "20230808T195901_20230808T195945_T10WEE",
    "20230811T200901_20230811T200857_T10WEE",
    "20230812T202849_20230812T202934_T10WEE",
    "20230813T195859_20230813T195900_T10WEE",
    "20230814T201851_20230814T202428_T10WEE",
    "20230809T201859_20230809T201853_T10WEE", # <- this one is broken
]

for i,poly in enumerate(bounds):
    print(f"Polygon {i}:")
    for image_id in images:
        try:
            ( geedim.MaskedImage
                .from_id("COPERNICUS/S2_SR_HARMONIZED/" + image_id, region=poly)
                .ee_image.getInfo()
            )
            print(f" {image_id} -> √")
        except ee.ee_exception.EEException as e:
            print(f" {image_id} -> error: {e}")
```

which produces

```
Polygon 0:
 20230801T200901_20230801T200857_T10WEE -> √
 20230802T202849_20230802T202847_T10WEE -> √
 20230803T195859_20230803T200054_T10WEE -> √
 20230804T201851_20230804T201853_T10WEE -> √
 20230804T201851_20230804T202447_T10WEE -> √
 20230806T200859_20230806T201056_T10WEE -> √
 20230807T202851_20230807T203151_T10WEE -> √
 20230808T195901_20230808T195945_T10WEE -> √
 20230811T200901_20230811T200857_T10WEE -> √
 20230812T202849_20230812T202934_T10WEE -> √
 20230813T195859_20230813T195900_T10WEE -> √
 20230814T201851_20230814T202428_T10WEE -> √
 20230809T201859_20230809T201853_T10WEE -> error: Image.rename: Parameter 'input' is required.
Polygon 1:
 20230801T200901_20230801T200857_T10WEE -> √
 20230802T202849_20230802T202847_T10WEE -> √
 20230803T195859_20230803T200054_T10WEE -> √
 20230804T201851_20230804T201853_T10WEE -> √
 20230804T201851_20230804T202447_T10WEE -> √
 20230806T200859_20230806T201056_T10WEE -> √
 20230807T202851_20230807T203151_T10WEE -> √
 20230808T195901_20230808T195945_T10WEE -> √
 20230811T200901_20230811T200857_T10WEE -> √
 20230812T202849_20230812T202934_T10WEE -> √
 20230813T195859_20230813T195900_T10WEE -> √
 20230814T201851_20230814T202428_T10WEE -> √
 20230809T201859_20230809T201853_T10WEE -> error: Image.rename: Parameter 'input' is required.
Polygon 2:
 20230801T200901_20230801T200857_T10WEE -> √
 20230802T202849_20230802T202847_T10WEE -> √
 20230803T195859_20230803T200054_T10WEE -> √
 20230804T201851_20230804T201853_T10WEE -> √
 20230804T201851_20230804T202447_T10WEE -> √
 20230806T200859_20230806T201056_T10WEE -> √
 20230807T202851_20230807T203151_T10WEE -> √
 20230808T195901_20230808T195945_T10WEE -> √
 20230811T200901_20230811T200857_T10WEE -> √
 20230812T202849_20230812T202934_T10WEE -> √
 20230813T195859_20230813T195900_T10WEE -> √
 20230814T201851_20230814T202428_T10WEE -> √
 20230809T201859_20230809T201853_T10WEE -> error: Image.rename: Parameter 'input' is required.
```





Discussion / hints:
Thanks for raising this.  The ``COPERNICUS/S2_CLOUD_PROBABILITY/20230809T201859_20230809T201853_T10WEE`` image geedim needs for cloud masking ``COPERNICUS/S2_SR_HARMONIZED/20230809T201859_20230809T201853_T10WEE`` doesn't exist.  I am looking into a way of dealing with this situation.

Fix patch:
diff --git a/geedim/download.py b/geedim/download.py
index c145e1b..4d68a59 100644
--- a/geedim/download.py
+++ b/geedim/download.py
@@ -52,7 +52,7 @@
 
 
 class BaseImage:
-    _float_nodata = float('nan')
+    _float_nodata = float('-inf')
     _desc_width = 50
     _default_resampling = ResamplingMethod.near
     _ee_max_tile_size = 32
@@ -984,6 +984,9 @@ def download(
 
             def download_tile(tile):
                 """Download a tile and write into the destination GeoTIFF."""
+                # Note: due to an Earth Engine bug (https://issuetracker.google.com/issues/350528377), tile_array is
+                # float64 for uint32 and int64 exp_image types.  The tile_array nodata value is as it should be
+                # though, so conversion to the exp_image / GeoTIFF type happens ok below.
                 tile_array = tile.download(session=session, bar=bar)
                 with out_lock:
                     out_ds.write(tile_array, window=tile.window)
diff --git a/geedim/mask.py b/geedim/mask.py
index f7bf5f7..20344e7 100644
--- a/geedim/mask.py
+++ b/geedim/mask.py
@@ -13,14 +13,16 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
+
 import logging
 from typing import Dict
 
 import ee
+
 import geedim.schema
 from geedim.download import BaseImage
 from geedim.enums import CloudMaskMethod
-from geedim.utils import split_id, get_projection
+from geedim.utils import get_projection, split_id
 
 logger = logging.getLogger(__name__)
 
@@ -71,15 +73,23 @@ def __init__(self, ee_image: ee.Image, mask: bool = _default_mask, region: dict
             Maximum distance (m) to look for clouds when forming the 'cloud distance' band.  Valid for
             Sentinel-2 images.
         """
-        # TODO: consider adding proj_scale parameter here, rather than in _set_region_stats, then it can be re-used in
-        #  S2 cloud masking and distance
         BaseImage.__init__(self, ee_image)
+        self._ee_projection = None
         self._add_aux_bands(**kwargs)  # add any mask and cloud distance bands
         if region:
             self._set_region_stats(region)
         if mask:
             self.mask_clouds()
 
+    @property
+    def _ee_proj(self) -> ee.Projection:
+        """Projection to use for mask calculations and statistics."""
+        if self._ee_projection is None:
+            # use the minimum scale projection for the base class / generic case (excludes non-fixed projections with
+            # 1 deg. scales)
+            self._ee_projection = get_projection(self._ee_image, min_scale=True)
+        return self._ee_projection
+
     @staticmethod
     def from_id(image_id: str, **kwargs) -> 'MaskedImage':
         """
@@ -97,7 +107,6 @@ def from_id(image_id: str, **kwargs) -> 'MaskedImage':
         MaskedImage
             A MaskedImage, or sub-class instance.
         """
-
         cls = class_from_id(image_id)
         ee_image = ee.Image(image_id)
         gd_image = cls(ee_image, **kwargs)
@@ -138,64 +147,61 @@ def _set_region_stats(self, region: Dict = None, scale: float = None):
         region : dict, ee.Geometry, optional
             Region inside of which to find statistics.  If not specified, the image footprint is used.
         scale: float, optional
-            Re-project to this scale when finding statistics.
+            Re-project to this scale when finding statistics.  Defaults to the scale of
+            :attr:`~MaskedImage._ee_proj`.  Should be provided if the encapsulated image is a composite / without a
+            fixed projection.
         """
         if not region:
             region = self.ee_image.geometry()  # use the image footprint
 
-        proj = get_projection(self.ee_image, min_scale=False)  # get projection of minimum scale band
-        # If _proj_scale is set, use that as the scale, otherwise use the proj.nomimalScale().  For non-composite images
-        # these should be the same value.  For composite images, there is no `fixed` projection, hence the
-        # need for _proj_scale.
-        scale = scale or proj.nominalScale()
+        # composite images have no fixed projection (scale = 1deg) and need the scale kwarg
+        scale = scale or self._ee_proj.nominalScale()
 
         # Find the fill portion as the (sum over the region of FILL_MASK) divided by (sum over the region of a constant
-        # image (==1)).  We take this approach rather than using a mean reducer, as this does not find the mean over
-        # the region, but the mean over the part of the region covered by the image.
-        stats_image = ee.Image(
-            [self.ee_image.select('FILL_MASK').unmask(), ee.Image(1).rename('REGION_SUM')]
-        )  # yapf: disable
+        # image (==1)).  Note that a mean reducer, does not find the mean over the region, but the mean over the part
+        # of the region covered by the image.
+        stats_image = ee.Image([self.ee_image.select('FILL_MASK').unmask(), ee.Image(1).rename('REGION_SUM')])
         # Note: sometimes proj has no EPSG in crs(), hence use crs=proj and not crs=proj.crs() below
         sums = stats_image.reduceRegion(
-            reducer="sum", geometry=region, crs=proj, scale=scale, bestEffort=True, maxPixels=1e6
+            reducer="sum", geometry=region, crs=self._ee_proj, scale=scale, bestEffort=True, maxPixels=1e6
         )
 
-        fill_portion = (
-            ee.Number(sums.get('FILL_MASK')).divide(ee.Number(sums.get('REGION_SUM'))).multiply(100)
-        )
+        fill_portion = ee.Number(sums.get('FILL_MASK')).divide(ee.Number(sums.get('REGION_SUM'))).multiply(100)
 
         # set the encapsulated image properties
         self.ee_image = self.ee_image.set('FILL_PORTION', fill_portion)
         # set CLOUDLESS_PORTION=100 for the generic case, where cloud/shadow masking is not supported
-        self.ee_image = self.ee_image.set('CLOUDLESS_PORTION', 100.)
+        self.ee_image = self.ee_image.set('CLOUDLESS_PORTION', 100.0)
 
     def mask_clouds(self):
-        """ Apply the cloud/shadow mask if supported, otherwise apply the fill mask. """
+        """Apply the cloud/shadow mask if supported, otherwise apply the fill mask."""
         self.ee_image = self.ee_image.updateMask(self.ee_image.select('FILL_MASK'))
 
 
 class CloudMaskedImage(MaskedImage):
-    """ A base class for encapsulating cloud/shadow masked images. """
+    """A base class for encapsulating cloud/shadow masked images."""
 
     def _cloud_dist(self, cloudless_mask: ee.Image = None, max_cloud_dist: float = 5000) -> ee.Image:
-        """ Find the cloud/shadow distance in units of 10m. """
+        """Find the cloud/shadow distance in units of 10m."""
         if not cloudless_mask:
             cloudless_mask = self.ee_image.select('CLOUDLESS_MASK')
-        proj = get_projection(self.ee_image, min_scale=False)  # use maximum scale projection to save processing time
 
-        # Note that initial *MASK bands before any call to mask_clouds(), are themselves masked, so this cloud/shadow
-        # mask excludes (i.e. masks) already masked pixels.  This avoids finding distance to e.g. scanline errors in
-        # Landsat-7.
         cloud_shadow_mask = cloudless_mask.Not()
-        cloud_pix = ee.Number(max_cloud_dist).divide(proj.nominalScale()).round()  # cloud_dist in pixels
-
-        # Find distance to nearest cloud/shadow (units of 10m).
-        cloud_dist = cloud_shadow_mask.fastDistanceTransform(
-            neighborhood=cloud_pix, units='pixels', metric='squared_euclidean'
-        ).sqrt().multiply(proj.nominalScale().divide(10))
+        cloud_pix = ee.Number(max_cloud_dist).divide(self._ee_proj.nominalScale()).round()  # cloud_dist in pixels
+
+        # Find distance to nearest cloud/shadow (units of 10m).  Distances are found for all pixels, including masked /
+        # invalid pixels, which are treated as 0 (non cloud/shadow).
+        cloud_dist = (
+            cloud_shadow_mask.fastDistanceTransform(neighborhood=cloud_pix, units='pixels', metric='squared_euclidean')
+            .sqrt()
+            .multiply(self._ee_proj.nominalScale().divide(10))
+        )
 
         # Reproject to force calculation at correct scale.
-        cloud_dist = cloud_dist.reproject(crs=proj, scale=proj.nominalScale()).rename('CLOUD_DIST')
+        cloud_dist = cloud_dist.reproject(crs=self._ee_proj, scale=self._ee_proj.nominalScale()).rename('CLOUD_DIST')
+
+        # Prevent use of invalid pixels.
+        cloud_dist = cloud_dist.updateMask(cloudless_mask.mask())
 
         # Clip cloud_dist to max_cloud_dist.
         cloud_dist = cloud_dist.where(cloud_dist.gt(ee.Image(max_cloud_dist / 10)), max_cloud_dist / 10)
@@ -213,30 +219,27 @@ def _set_region_stats(self, region: Dict = None, scale: float = None):
         region : dict, ee.Geometry, optional
             Region inside of which to find statistics.  If not specified, the image footprint is used.
         scale: float, optional
-            Re-project to this scale when finding statistics.
+            Re-project to this scale when finding statistics.  Defaults to the scale of
+            :attr:`~MaskedImage._ee_proj`.  Should be provided if the encapsulated image is a composite / without a
+            fixed projection.
         """
         if not region:
             region = self.ee_image.geometry()  # use the image footprint
 
-        proj = get_projection(self.ee_image, min_scale=False)  # get projection of minimum scale band
-        # If _proj_scale is set, use that as the scale, otherwise use the proj.nomimalScale().  For non-composite images
-        # these should be the same value.  For composite images, there is no `fixed` projection, hence the
-        # need for _proj_scale.
-        scale = scale or proj.nominalScale()
+        # composite images have no fixed projection (scale = 1deg) and need the scale kwarg
+        scale = scale or self._ee_proj.nominalScale()
 
-        # Find the fill portion as the (sum over the region of FILL_MASK) divided by (sum over the region of a constant
-        # image (==1)).  We take this approach rather than using a mean reducer, as this does not find the mean over
-        # the region, but the mean over the part of the region covered by the image.  Then we find cloudless portion
-        # as portion of fill that is cloudless.
+        # Find the fill portion as the (sum over the region of FILL_MASK) divided by (sum over the region of a
+        # constant image (==1)).  Then find cloudless portion as portion of fill that is cloudless.  Note that a mean
+        # reducer, does not find the mean over the region, but the mean over the part of the region covered by the
+        # image.
         stats_image = ee.Image(
             [self.ee_image.select(['FILL_MASK', 'CLOUDLESS_MASK']).unmask(), ee.Image(1).rename('REGION_SUM')]
-        )  # yapf: disable
+        )
         sums = stats_image.reduceRegion(
-            reducer="sum", geometry=region, crs=proj, scale=scale, bestEffort=True, maxPixels=1e6
+            reducer="sum", geometry=region, crs=self._ee_proj, scale=scale, bestEffort=True, maxPixels=1e6
         ).rename(['FILL_MASK', 'CLOUDLESS_MASK'], ['FILL_PORTION', 'CLOUDLESS_PORTION'])
-        fill_portion = (
-            ee.Number(sums.get('FILL_PORTION')).divide(ee.Number(sums.get('REGION_SUM'))).multiply(100)
-        )
+        fill_portion = ee.Number(sums.get('FILL_PORTION')).divide(ee.Number(sums.get('REGION_SUM'))).multiply(100)
         cloudless_portion = (
             ee.Number(sums.get('CLOUDLESS_PORTION')).divide(ee.Number(sums.get('FILL_PORTION'))).multiply(100)
         )
@@ -255,7 +258,7 @@ def _aux_image(self, **kwargs) -> ee.Image:
         raise NotImplementedError('This virtual method should be overridden by derived classes.')
 
     def mask_clouds(self):
-        """ Apply the cloud/shadow mask. """
+        """Apply the cloud/shadow mask."""
         self.ee_image = self.ee_image.updateMask(self.ee_image.select('CLOUDLESS_MASK'))
 
 
@@ -271,6 +274,13 @@ class LandsatImage(CloudMaskedImage):
     * LANDSAT/LC09/C02/T1_L2
     """
 
+    @property
+    def _ee_proj(self) -> ee.Projection:
+        if self._ee_projection is None:
+            # use the default projection (all landsat bands have same scale)
+            self._ee_projection = self._ee_image.projection()
+        return self._ee_projection
+
     def _aux_image(self, mask_shadows: bool = True, mask_cirrus: bool = True, max_cloud_dist: int = 5000) -> ee.Image:
         """
         Retrieve the auxiliary image containing cloud/shadow masks and cloud distance.
@@ -315,12 +325,28 @@ def _aux_image(self, mask_shadows: bool = True, mask_cirrus: bool = True, max_cl
 
 
 class Sentinel2ClImage(CloudMaskedImage):
-    """ Base class for cloud/shadow masking of Sentinel-2 TOA and SR images. """
+    """Base class for cloud/shadow masking of Sentinel-2 TOA and SR images."""
+
+    @property
+    def _ee_proj(self) -> ee.Projection:
+        if self._ee_projection is None:
+            # use the B1 projection with maximum scale (60m) to reduce processing times (some S2 images have empty QA
+            # bands with no fixed projection, so utils.get_projection(min_scale=False) should not be used here).
+            self._ee_projection = self._ee_image.select(0).projection()
+        return self._ee_projection
 
     def _aux_image(
-        self, s2_toa: bool = False, mask_cirrus: bool = True, mask_shadows: bool = True,
-        mask_method: CloudMaskMethod = CloudMaskMethod.cloud_prob, prob: float = 60, dark: float = 0.15,
-        shadow_dist: int = 1000, buffer: int = 50, cdi_thresh: float = None, max_cloud_dist: int = 5000
+        self,
+        s2_toa: bool = False,
+        mask_cirrus: bool = True,
+        mask_shadows: bool = True,
+        mask_method: CloudMaskMethod = CloudMaskMethod.cloud_prob,
+        prob: float = 60,
+        dark: float = 0.15,
+        shadow_dist: int = 1000,
+        buffer: int = 50,
+        cdi_thresh: float = None,
+        max_cloud_dist: int = 5000,
     ) -> ee.Image:
         """
         Derive cloud, shadow and validity masks for the encapsulated image.
@@ -362,12 +388,23 @@ def _aux_image(
         ee.Image
             An Earth Engine image containing *_MASK and CLOUD_DIST bands.
         """
+        # TODO: add warning about post Feb 2022 validity when qa method is used
         mask_method = CloudMaskMethod(mask_method)
 
         def get_cloud_prob(ee_im):
-            """Get the cloud probability image from COPERNICUS/S2_CLOUD_PROBABILITY that corresponds to `ee_im`."""
+            """Get the cloud probability image from COPERNICUS/S2_CLOUD_PROBABILITY that corresponds to `ee_im`, if it
+            exists.  Otherwise, get a 100% probability.
+            """
+            # default masked cloud probability image
+            default_cloud_prob = ee.Image().updateMask(0)
+
+            # get the cloud probability image if it exists, otherwise revert to default_cloud_prob
+            # (ee.Image.linkCollection() has a similar functionality but the masked image it returns when cloud
+            # probability doesn't exist is 0 rather than nodata on download)
             filt = ee.Filter.eq('system:index', ee_im.get('system:index'))
             cloud_prob = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY').filter(filt).first()
+            cloud_prob = ee.Image(ee.List([cloud_prob, default_cloud_prob]).reduce(ee.Reducer.firstNonNull()))
+
             return cloud_prob.rename('CLOUD_PROB')
 
         def get_cloud_mask(ee_im, cloud_prob=None):
@@ -377,10 +414,15 @@ def get_cloud_mask(ee_im, cloud_prob=None):
                     cloud_prob = get_cloud_prob(ee_im)
                 cloud_mask = cloud_prob.gte(prob)
             else:
-                qa = ee_im.select('QA60')
+                # mask QA60 if it is post Feb 2022 to invalidate the cloud mask (post Feb 2022 QA60 is no longer
+                # populated and may be masked / transparent or zero).
+                qa_valid = ee.Number(ee_im.date().difference(ee.Date('2022-02-01'), 'days').lt(0))
+                qa = ee_im.select('QA60').updateMask(qa_valid)
+
                 cloud_mask = qa.bitwiseAnd(1 << 10).neq(0)
                 if mask_cirrus:
                     cloud_mask = cloud_mask.Or(qa.bitwiseAnd(1 << 11).neq(0))
+
             return cloud_mask.rename('CLOUD_MASK')
 
         def get_cdi_cloud_mask(ee_im):
@@ -393,7 +435,7 @@ def get_cdi_cloud_mask(ee_im):
             else:
                 # get the Sentinel-2 TOA image that corresponds to ee_im
                 idx = ee_im.get('system:index')
-                s2_toa_image = (ee.ImageCollection('COPERNICUS/S2').filter(ee.Filter.eq('system:index', idx)).first())
+                s2_toa_image = ee.ImageCollection('COPERNICUS/S2').filter(ee.Filter.eq('system:index', idx)).first()
             cdi_image = ee.Algorithms.Sentinel2.CDI(s2_toa_image)
             return cdi_image.lt(cdi_thresh).rename('CDI_CLOUD_MASK')
 
@@ -401,29 +443,39 @@ def get_shadow_mask(ee_im, cloud_mask):
             """Given a cloud mask, get a shadow mask for ee_im."""
             dark_mask = ee_im.select('B8').lt(dark * 1e4)
             if not s2_toa:
+                # exclude water
                 dark_mask = ee_im.select('SCL').neq(6).And(dark_mask)
 
-            proj = get_projection(ee_im, min_scale=False)
             # Note:
             # S2 MEAN_SOLAR_AZIMUTH_ANGLE (SAA) appears to be measured clockwise with 0 at N (i.e. shadow goes in the
             # opposite direction), directionalDistanceTransform() angle appears to be measured clockwise with 0 at W.
             # So we need to add/subtract 180 to SAA to get shadow angle in S2 convention, then add 90 to get
-            # directionalDistanceTransform() convention i.e. we need to add -180 + 90 = -90 to the SAA.  This is not
+            # directionalDistanceTransform() convention i.e. we need to add -180 + 90 = -90 to SAA.  This is not
             # the same as in the EE tutorial which is 90-SAA
             # (https://developers.google.com/earth-engine/tutorials/community/sentinel-2-s2cloudless).
             shadow_azimuth = ee.Number(-90).add(ee.Number(ee_im.get('MEAN_SOLAR_AZIMUTH_ANGLE')))
-            proj_pixels = ee.Number(shadow_dist).divide(proj.nominalScale()).round()
-            # Project the cloud mask in the direction of the shadows it will cast.
+            proj_pixels = ee.Number(shadow_dist).divide(self._ee_proj.nominalScale()).round()
+
+            # Project the cloud mask in the direction of the shadows it will cast (can project the mask into invalid
+            # areas).
             cloud_cast_proj = cloud_mask.directionalDistanceTransform(shadow_azimuth, proj_pixels).select('distance')
-            # The reproject is necessary to force calculation at the correct scale - the coarse proj scale is used to
-            # improve processing times.
-            cloud_cast_mask = cloud_cast_proj.mask().reproject(crs=proj.crs(), scale=proj.nominalScale())
-            return cloud_cast_mask.And(dark_mask).rename('SHADOW_MASK')
 
-        # gather and combine the various masks
+            # Reproject to force calculation at the correct scale.
+            cloud_cast_mask = cloud_cast_proj.mask().reproject(crs=self._ee_proj, scale=self._ee_proj.nominalScale())
+
+            # Remove any projections in invalid areas.
+            cloud_cast_mask = cloud_cast_mask.updateMask(cloud_mask.mask())
+
+            # Find shadow mask as intersection between projected clouds and dark areas.
+            shadow_mask = cloud_cast_mask.And(dark_mask)
+
+            return shadow_mask.rename('SHADOW_MASK')
+
+        # combine the various masks
         ee_image = self.ee_image
         cloud_prob = get_cloud_prob(ee_image) if mask_method == CloudMaskMethod.cloud_prob else None
         cloud_mask = get_cloud_mask(ee_image, cloud_prob=cloud_prob)
+
         if cdi_thresh is not None:
             cloud_mask = cloud_mask.And(get_cdi_cloud_mask(ee_image))
         if mask_shadows:
@@ -434,9 +486,11 @@ def get_shadow_mask(ee_im, cloud_mask):
 
         # do a morphological opening type operation that removes small (20m) blobs from the mask and then dilates
         cloud_shadow_mask = cloud_shadow_mask.focal_min(20, units='meters').focal_max(buffer, units='meters')
+
         # derive a fill mask from the Earth Engine mask for the surface reflectance bands
         fill_mask = ee_image.select('B.*').mask().reduce(ee.Reducer.allNonZero())
-        # Clip this mask to the image footprint.  (Without this step we get memory limit errors on download.)
+
+        # clip this mask to the image footprint (without this step we get memory limit errors on download)
         fill_mask = fill_mask.clip(ee_image.geometry()).rename('FILL_MASK')
 
         # combine all masks into cloudless_mask
@@ -454,21 +508,21 @@ def get_shadow_mask(ee_im, cloud_mask):
 
 
 class Sentinel2SrClImage(Sentinel2ClImage):
-    """ Class for cloud/shadow masking of Sentinel-2 SR (COPERNICUS/S2_SR) images. """
+    """Class for cloud/shadow masking of Sentinel-2 SR (COPERNICUS/S2_SR) images."""
 
     def _aux_image(self, s2_toa: bool = False, **kwargs) -> ee.Image:
         return Sentinel2ClImage._aux_image(self, s2_toa=False, **kwargs)
 
 
 class Sentinel2ToaClImage(Sentinel2ClImage):
-    """ Class for cloud/shadow masking of Sentinel-2 TOA (COPERNICUS/S2) images. """
+    """Class for cloud/shadow masking of Sentinel-2 TOA (COPERNICUS/S2) images."""
 
     def _aux_image(self, s2_toa: bool = False, **kwargs) -> ee.Image:
         return Sentinel2ClImage._aux_image(self, s2_toa=True, **kwargs)
 
 
 def class_from_id(image_id: str) -> type:
-    """ Return the *Image class that corresponds to the provided Earth Engine image/collection ID. """
+    """Return the *Image class that corresponds to the provided Earth Engine image/collection ID."""
     ee_coll_name, _ = split_id(image_id)
     if image_id in geedim.schema.collection_schema:
         return geedim.schema.collection_schema[image_id]['image_type']
diff --git a/geedim/tile.py b/geedim/tile.py
index d683ad1..39e98d8 100644
--- a/geedim/tile.py
+++ b/geedim/tile.py
@@ -127,11 +127,6 @@ def _download_to_array(self, url: str, session: requests.Session = None, bar: tq
         with utils.suppress_rio_logs(), env, MemoryFile(gtiff_bytes) as mem_file:
             with mem_file.open() as ds:
                 array = ds.read()
-                if (array.dtype == np.dtype('float32')) or (array.dtype == np.dtype('float64')):
-                    # GEE sets nodata to -inf for float data types, (but does not populate the nodata field).
-                    # rasterio won't allow nodata=-inf, so this is a workaround to change nodata to nan at
-                    # source.
-                    array[np.isinf(array)] = np.nan
         return array
 
     def download(
diff --git a/geedim/utils.py b/geedim/utils.py
index c5f00c9..275526f 100644
--- a/geedim/utils.py
+++ b/geedim/utils.py
@@ -23,7 +23,7 @@
 import time
 from contextlib import contextmanager
 from threading import Thread
-from typing import Tuple, Optional
+from typing import Optional, Tuple
 
 import ee
 import numpy as np
@@ -164,7 +164,7 @@ def get_bounds(filename: pathlib.Path, expand: float = 5):
     return src_bbox_wgs84
 
 
-def get_projection(image, min_scale=True):
+def get_projection(image: ee.Image, min_scale: bool = True) -> ee.Projection:
     """
     Get the min/max scale projection of image bands.  Server side - no calls to getInfo().
     Adapted from from https://github.com/gee-community/gee_tools, MIT license.
@@ -186,23 +186,11 @@ def get_projection(image, min_scale=True):
         raise TypeError('image is not an instance of ee.Image')
 
     bands = image.bandNames()
+    scales = bands.map(lambda band: image.select(ee.String(band)).projection().nominalScale())
+    projs = bands.map(lambda band: image.select(ee.String(band)).projection())
+    projs = projs.sort(scales)
 
-    compare = ee.Number.lte if min_scale else ee.Number.gte
-    init_proj = image.select(0).projection()
-
-    def compare_scale(name, prev_proj):
-        """Server side comparison of band scales"""
-        prev_proj = ee.Projection(prev_proj)
-        prev_scale = prev_proj.nominalScale()
-
-        curr_proj = image.select([name]).projection()
-        curr_scale = ee.Number(curr_proj.nominalScale())
-
-        condition = compare(curr_scale, prev_scale)
-        comp_proj = ee.Algorithms.If(condition, curr_proj, prev_proj)
-        return ee.Projection(comp_proj)
-
-    return ee.Projection(bands.iterate(compare_scale, init_proj))
+    return ee.Projection(projs.get(0) if min_scale else projs.get(-1))
 
 
 class Spinner(Thread):
diff --git a/pyproject.toml b/pyproject.toml
index aa48457..2344c30 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -5,7 +5,7 @@ readme = 'README.rst'
 requires-python = '>=3.6'
 dependencies = [
     'numpy>=1.19',
-    'rasterio>=1.1',
+    'rasterio>=1.3',
     'click>=8',
     'tqdm>=4.6',
     'earthengine-api>=0.1.2',


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 6d23ed4..887034e 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -13,13 +13,16 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
+
 import pathlib
 from typing import Dict, List
 
 import ee
 import pytest
 from click.testing import CliRunner
+
 from geedim import Initialize, MaskedImage
+from geedim.mask import Sentinel2SrClImage
 from geedim.utils import root_path
 
 
@@ -31,31 +34,34 @@ def ee_init():
 
 @pytest.fixture(scope='session')
 def region_25ha() -> Dict:
-    """ A geojson polygon defining a 500x500m region. """
+    """A geojson polygon defining a 500x500m region."""
     return {
         "type": "Polygon",
-        "coordinates":
-            [[[21.6389, -33.4520], [21.6389, -33.4474], [21.6442, -33.4474], [21.6442, -33.4520], [21.6389, -33.4520]]]
+        "coordinates": [
+            [[21.6389, -33.4520], [21.6389, -33.4474], [21.6442, -33.4474], [21.6442, -33.4520], [21.6389, -33.4520]]
+        ],
     }
 
 
 @pytest.fixture(scope='session')
 def region_100ha() -> Dict:
-    """ A geojson polygon defining a 1x1km region. """
+    """A geojson polygon defining a 1x1km region."""
     return {
         "type": "Polygon",
-        "coordinates":
-            [[[21.6374, -33.4547], [21.6374, -33.4455], [21.6480, -33.4455], [21.6480, -33.4547], [21.6374, -33.4547]]]
+        "coordinates": [
+            [[21.6374, -33.4547], [21.6374, -33.4455], [21.6480, -33.4455], [21.6480, -33.4547], [21.6374, -33.4547]]
+        ],
     }
 
 
 @pytest.fixture(scope='session')
 def region_10000ha() -> Dict:
-    """ A geojson polygon defining a 10x10km region. """
+    """A geojson polygon defining a 10x10km region."""
     return {
         "type": "Polygon",
-        "coordinates":
-            [[[21.5893, -33.4964], [21.5893, -33.4038], [21.6960, -33.4038], [21.6960, -33.4964], [21.5893, -33.4964]]]
+        "coordinates": [
+            [[21.5893, -33.4964], [21.5893, -33.4038], [21.6960, -33.4038], [21.6960, -33.4964], [21.5893, -33.4964]]
+        ],
     }
 
 
@@ -66,75 +72,91 @@ def const_image_25ha_file() -> pathlib.Path:
 
 @pytest.fixture(scope='session')
 def l4_image_id() -> str:
-    """ Landsat-4 EE ID for image that covers `region_*ha`, with partial cloud cover only for `region10000ha`.  """
+    """Landsat-4 EE ID for image that covers `region_*ha`, with partial cloud cover only for `region10000ha`."""
     return 'LANDSAT/LT04/C02/T1_L2/LT04_173083_19880310'
 
 
 @pytest.fixture(scope='session')
 def l5_image_id() -> str:
-    """ Landsat-5 EE ID for image that covers `region_*ha` with partial cloud cover.  """
+    """Landsat-5 EE ID for image that covers `region_*ha` with partial cloud cover."""
     return 'LANDSAT/LT05/C02/T1_L2/LT05_173083_20051112'  # 'LANDSAT/LT05/C02/T1_L2/LT05_173083_20070307'
 
 
 @pytest.fixture(scope='session')
 def l7_image_id() -> str:
-    """ Landsat-7 EE ID for image that covers `region_*ha` with partial cloud cover.  """
+    """Landsat-7 EE ID for image that covers `region_*ha` with partial cloud cover."""
     return 'LANDSAT/LE07/C02/T1_L2/LE07_173083_20220119'  # 'LANDSAT/LE07/C02/T1_L2/LE07_173083_20200521'
 
 
 @pytest.fixture(scope='session')
 def l8_image_id() -> str:
-    """ Landsat-8 EE ID for image that covers `region_*ha` with partial cloud cover.  """
+    """Landsat-8 EE ID for image that covers `region_*ha` with partial cloud cover."""
     return 'LANDSAT/LC08/C02/T1_L2/LC08_173083_20180217'  # 'LANDSAT/LC08/C02/T1_L2/LC08_173083_20171113'
 
 
 @pytest.fixture(scope='session')
 def l9_image_id() -> str:
-    """ Landsat-9 EE ID for image that covers `region_*ha` with partial cloud cover. """
+    """Landsat-9 EE ID for image that covers `region_*ha` with partial cloud cover."""
     return 'LANDSAT/LC09/C02/T1_L2/LC09_173083_20220308'  # 'LANDSAT/LC09/C02/T1_L2/LC09_173083_20220103'
 
 
 @pytest.fixture(scope='session')
 def landsat_image_ids(l4_image_id, l5_image_id, l7_image_id, l8_image_id, l9_image_id) -> List[str]:
-    """ Landsat4-9 EE IDs for images that covers `region_*ha` with partial cloud cover. """
+    """Landsat4-9 EE IDs for images that covers `region_*ha` with partial cloud cover."""
     return [l4_image_id, l5_image_id, l7_image_id, l8_image_id, l9_image_id]
 
 
 @pytest.fixture(scope='session')
 def s2_sr_image_id() -> str:
-    """ Sentinel-2 SR EE ID for image that covers `region_*ha` with partial cloud cover. """
+    """Sentinel-2 SR EE ID for image that covers `region_*ha` with partial cloud cover."""
     # 'COPERNICUS/S2/20220107T081229_20220107T083059_T34HEJ'
     return 'COPERNICUS/S2_SR/20211004T080801_20211004T083709_T34HEJ'
 
 
 @pytest.fixture(scope='session')
 def s2_toa_image_id() -> str:
-    """ Sentinel-2 TOA EE ID for image that covers `region_*ha` with partial cloud cover. """
+    """Sentinel-2 TOA EE ID for image that covers `region_*ha` with partial cloud cover."""
     return 'COPERNICUS/S2/20220107T081229_20220107T083059_T34HEJ'
 
 
 @pytest.fixture(scope='session')
 def s2_sr_hm_image_id() -> str:
-    """ Harmonised Sentinel-2 SR EE ID for image that covers `region_*ha` with partial cloud cover. """
+    """Harmonised Sentinel-2 SR EE ID for image that covers `region_*ha` with partial cloud cover."""
     # 'COPERNICUS/S2/20220107T081229_20220107T083059_T34HEJ'
     return 'COPERNICUS/S2_SR_HARMONIZED/20211004T080801_20211004T083709_T34HEJ'
 
 
+@pytest.fixture(scope='session')
+def s2_sr_hm_qa_mask_image_id() -> str:
+    """Harmonised Sentinel-2 SR EE ID for image with masked QA* data, that covers `region_*ha` with partial cloud
+    cover.
+    """
+    return 'COPERNICUS/S2_SR_HARMONIZED/20240426T080609_20240426T083054_T34HEJ'
+
+
+@pytest.fixture(scope='session')
+def s2_sr_hm_qa_zero_image_id() -> str:
+    """Harmonised Sentinel-2 SR EE ID for image with zero QA* data, that covers `region_*ha` with partial cloud
+    cover.
+    """
+    return 'COPERNICUS/S2_SR_HARMONIZED/20230407T080611_20230407T083613_T34HEJ'
+
+
 @pytest.fixture(scope='session')
 def s2_toa_hm_image_id() -> str:
-    """ Harmonised Sentinel-2 TOA EE ID for image that covers `region_*ha` with partial cloud cover. """
+    """Harmonised Sentinel-2 TOA EE ID for image that covers `region_*ha` with partial cloud cover."""
     return 'COPERNICUS/S2_HARMONIZED/20220107T081229_20220107T083059_T34HEJ'
 
 
 @pytest.fixture(scope='session')
 def s2_image_ids(s2_sr_image_id, s2_toa_image_id, s2_sr_hm_image_id, s2_toa_hm_image_id) -> List[str]:
-    """ Sentinel-2 TOA/SR EE IDs for images that covers `region_*ha` with partial cloud cover. """
+    """Sentinel-2 TOA/SR EE IDs for images that covers `region_*ha` with partial cloud cover."""
     return [s2_sr_image_id, s2_toa_image_id, s2_sr_hm_image_id, s2_toa_hm_image_id]
 
 
 @pytest.fixture(scope='session')
 def modis_nbar_image_id() -> str:
-    """ Global MODIS NBAR image ID.  WGS84 @ 500m. """
+    """Global MODIS NBAR image ID.  WGS84 @ 500m."""
     return 'MODIS/006/MCD43A4/2022_01_01'
 
 
@@ -149,25 +171,25 @@ def gch_image_id() -> str:
 
 @pytest.fixture(scope='session')
 def s1_sar_image_id() -> str:
-    """ Sentinel-1 SAR GRD EE image ID.  10m. """
+    """Sentinel-1 SAR GRD EE image ID.  10m."""
     return 'COPERNICUS/S1_GRD/S1A_IW_GRDH_1SDV_20220112T171750_20220112T171815_041430_04ED28_0A04'
 
 
 @pytest.fixture(scope='session')
 def gedi_agb_image_id() -> str:
-    """ GEDI aboveground biomass density EE image ID.  1km."""
+    """GEDI aboveground biomass density EE image ID.  1km."""
     return 'LARSE/GEDI/GEDI04_B_002'
 
 
 @pytest.fixture(scope='session')
 def gedi_cth_image_id() -> str:
-    """ GEDI canopy top height EE image ID.  25m."""
+    """GEDI canopy top height EE image ID.  25m."""
     return 'LARSE/GEDI/GEDI02_A_002_MONTHLY/202009_018E_036S'
 
 
 @pytest.fixture(scope='session')
 def landsat_ndvi_image_id() -> str:
-    """ Landsat 8-day NDVI composite EE image iD.  Composite in WGS84 with underlying 30m scale."""
+    """Landsat 8-day NDVI composite EE image iD.  Composite in WGS84 with underlying 30m scale."""
     return 'LANDSAT/LC08/C01/T1_8DAY_NDVI/20211219'
 
 
@@ -175,157 +197,193 @@ def landsat_ndvi_image_id() -> str:
 def generic_image_ids(
     modis_nbar_image_id, gch_image_id, s1_sar_image_id, gedi_agb_image_id, gedi_cth_image_id, landsat_ndvi_image_id
 ) -> List[str]:
-    """ A list of various EE image IDs for non-cloud/shadow masked images. """
+    """A list of various EE image IDs for non-cloud/shadow masked images."""
     return [
-        modis_nbar_image_id, gch_image_id, s1_sar_image_id, gedi_agb_image_id, gedi_cth_image_id, landsat_ndvi_image_id
+        modis_nbar_image_id,
+        gch_image_id,
+        s1_sar_image_id,
+        gedi_agb_image_id,
+        gedi_cth_image_id,
+        landsat_ndvi_image_id,
     ]
 
 
 @pytest.fixture(scope='session')
 def l4_masked_image(l4_image_id) -> MaskedImage:
-    """ Landsat-4 MaskedImage that covers `region_*ha`, with partial cloud cover only for `region10000ha`. """
+    """Landsat-4 MaskedImage that covers `region_*ha`, with partial cloud cover only for `region10000ha`."""
     return MaskedImage.from_id(l4_image_id)
 
 
 @pytest.fixture(scope='session')
 def l5_masked_image(l5_image_id) -> MaskedImage:
-    """ Landsat-5 MaskedImage that covers `region_*ha` with partial cloud cover. """
+    """Landsat-5 MaskedImage that covers `region_*ha` with partial cloud cover."""
     return MaskedImage.from_id(l5_image_id)
 
 
 @pytest.fixture(scope='session')
 def l7_masked_image(l7_image_id) -> MaskedImage:
-    """ Landsat-7 MaskedImage that covers `region_*ha` with partial cloud cover. """
+    """Landsat-7 MaskedImage that covers `region_*ha` with partial cloud cover."""
     return MaskedImage.from_id(l7_image_id)
 
 
 @pytest.fixture(scope='session')
 def l8_masked_image(l8_image_id) -> MaskedImage:
-    """ Landsat-8 MaskedImage that cover `region_*ha` with partial cloud cover. """
+    """Landsat-8 MaskedImage that cover `region_*ha` with partial cloud cover."""
     return MaskedImage.from_id(l8_image_id)
 
 
 @pytest.fixture(scope='session')
 def l9_masked_image(l9_image_id) -> MaskedImage:
-    """ Landsat-9 MaskedImage that covers `region_*ha` with partial cloud cover. """
+    """Landsat-9 MaskedImage that covers `region_*ha` with partial cloud cover."""
     return MaskedImage.from_id(l9_image_id)
 
 
 @pytest.fixture(scope='session')
 def landsat_masked_images(
-    l4_masked_image, l5_masked_image, l7_masked_image, l8_masked_image,
-    l9_masked_image
-) -> List[MaskedImage]: # yapf: disable
-    """ Landsat4-9 MaskedImage's that cover `region_*ha` with partial cloud cover. """
+    l4_masked_image, l5_masked_image, l7_masked_image, l8_masked_image, l9_masked_image
+) -> List[MaskedImage]:  # yapf: disable
+    """Landsat4-9 MaskedImage's that cover `region_*ha` with partial cloud cover."""
     return [l4_masked_image, l5_masked_image, l7_masked_image, l8_masked_image, l9_masked_image]
 
 
 @pytest.fixture(scope='session')
 def s2_sr_masked_image(s2_sr_image_id) -> MaskedImage:
-    """ Sentinel-2 SR MaskedImage that covers `region_*ha` with partial cloud cover. """
+    """Sentinel-2 SR MaskedImage that covers `region_*ha` with partial cloud cover."""
     return MaskedImage.from_id(s2_sr_image_id)
 
 
 @pytest.fixture(scope='session')
 def s2_toa_masked_image(s2_toa_image_id) -> MaskedImage:
-    """ Sentinel-2 TOA MaskedImage that covers `region_*ha` with partial cloud cover. """
+    """Sentinel-2 TOA MaskedImage that covers `region_*ha` with partial cloud cover."""
     return MaskedImage.from_id(s2_toa_image_id)
 
 
 @pytest.fixture(scope='session')
 def s2_sr_hm_masked_image(s2_sr_hm_image_id) -> MaskedImage:
-    """ Harmonised Sentinel-2 SR MaskedImage that covers `region_*ha` with partial cloud cover. """
+    """Harmonised Sentinel-2 SR MaskedImage that covers `region_*ha` with partial cloud cover."""
     return MaskedImage.from_id(s2_sr_hm_image_id)
 
 
+@pytest.fixture(scope='session')
+def s2_sr_hm_nocp_masked_image(s2_sr_hm_image_id) -> MaskedImage:
+    """Harmonised Sentinel-2 SR MaskedImage with no corresponding cloud probability, that covers `region_*ha` with
+    partial cloud cover.
+    """
+    # create an image with unknown id to prevent linking to cloud probability
+    ee_image = ee.Image(s2_sr_hm_image_id)
+    ee_image = ee_image.set('system:index', 'COPERNICUS/S2_HARMONIZED/unknown')
+    return Sentinel2SrClImage(ee_image, mask_method='cloud-prob')
+
+
+@pytest.fixture(scope='session')
+def s2_sr_hm_qa_mask_masked_image(s2_sr_hm_qa_mask_image_id: str) -> MaskedImage:
+    """Harmonised Sentinel-2 SR MaskedImage with masked QA* bands, that covers `region_*ha` with partial cloud cover."""
+    return MaskedImage.from_id(s2_sr_hm_qa_mask_image_id, mask_method='qa')
+
+
+@pytest.fixture(scope='session')
+def s2_sr_hm_qa_zero_masked_image(s2_sr_hm_qa_zero_image_id: str) -> MaskedImage:
+    """Harmonised Sentinel-2 SR MaskedImage with zero QA* bands, that covers `region_*ha` with partial cloud cover."""
+    return MaskedImage.from_id(s2_sr_hm_qa_zero_image_id, mask_method='qa')
+
+
 @pytest.fixture(scope='session')
 def s2_toa_hm_masked_image(s2_toa_hm_image_id) -> MaskedImage:
-    """ Harmonised Sentinel-2 TOA MaskedImage that covers `region_*ha` with partial cloud cover. """
+    """Harmonised Sentinel-2 TOA MaskedImage that covers `region_*ha` with partial cloud cover."""
     return MaskedImage.from_id(s2_toa_hm_image_id)
 
 
 @pytest.fixture(scope='session')
 def s2_masked_images(
     s2_sr_masked_image, s2_toa_masked_image, s2_sr_hm_masked_image, s2_toa_hm_masked_image
-) -> List[MaskedImage]: # yapf: disable
-    """ Sentinel-2 TOA and SR MaskedImage's that cover `region_*ha` with partial cloud cover. """
+) -> List[MaskedImage]:  # yapf: disable
+    """Sentinel-2 TOA and SR MaskedImage's that cover `region_*ha` with partial cloud cover."""
     return [s2_sr_masked_image, s2_toa_masked_image, s2_sr_hm_masked_image, s2_toa_hm_masked_image]
 
 
 @pytest.fixture(scope='session')
 def user_masked_image() -> MaskedImage:
-    """ A MaskedImage instance where the encapsulated image has no fixed projection or ID.  """
+    """A MaskedImage instance where the encapsulated image has no fixed projection or ID."""
     return MaskedImage(ee.Image([1, 2, 3]))
 
 
 @pytest.fixture(scope='session')
 def modis_nbar_masked_image(modis_nbar_image_id) -> MaskedImage:
-    """ MODIS NBAR MaskedImage with global coverage.  """
+    """MODIS NBAR MaskedImage with global coverage."""
     return MaskedImage.from_id(modis_nbar_image_id)
 
 
 @pytest.fixture(scope='session')
 def gch_masked_image(gch_image_id) -> MaskedImage:
-    """ Global Canopy Height (10m) MaskedImage. """
+    """Global Canopy Height (10m) MaskedImage."""
     return MaskedImage.from_id(gch_image_id)
 
 
 @pytest.fixture(scope='session')
 def s1_sar_masked_image(s1_sar_image_id) -> MaskedImage:
-    """ Sentinel-1 SAR GRD MaskedImage.  10m. """
+    """Sentinel-1 SAR GRD MaskedImage.  10m."""
     return MaskedImage.from_id(s1_sar_image_id)
 
 
 @pytest.fixture(scope='session')
 def gedi_agb_masked_image(gedi_agb_image_id) -> MaskedImage:
-    """ GEDI aboveground biomass density MaskedImage.  1km."""
+    """GEDI aboveground biomass density MaskedImage.  1km."""
     return MaskedImage.from_id(gedi_agb_image_id)
 
 
 @pytest.fixture(scope='session')
 def gedi_cth_masked_image(gedi_cth_image_id) -> MaskedImage:
-    """ GEDI canopy top height MaskedImage.  25m."""
+    """GEDI canopy top height MaskedImage.  25m."""
     return MaskedImage.from_id(gedi_cth_image_id)
 
 
 @pytest.fixture(scope='session')
 def landsat_ndvi_masked_image(landsat_ndvi_image_id) -> MaskedImage:
-    """ Landsat 8-day NDVI composite MaskedImage.  Composite in WGS84 with underlying 30m scale."""
+    """Landsat 8-day NDVI composite MaskedImage.  Composite in WGS84 with underlying 30m scale."""
     return MaskedImage.from_id(landsat_ndvi_image_id)
 
 
 @pytest.fixture(scope='session')
 def generic_masked_images(
-    modis_nbar_masked_image, gch_masked_image, s1_sar_masked_image, gedi_agb_masked_image, gedi_cth_masked_image,
-    landsat_ndvi_masked_image
+    modis_nbar_masked_image,
+    gch_masked_image,
+    s1_sar_masked_image,
+    gedi_agb_masked_image,
+    gedi_cth_masked_image,
+    landsat_ndvi_masked_image,
 ) -> List[MaskedImage]:
-    """ A list of various non-cloud/shadow MaskedImage's. """
+    """A list of various non-cloud/shadow MaskedImage's."""
     return [
-        modis_nbar_masked_image, gch_masked_image, s1_sar_masked_image, gedi_agb_masked_image, gedi_cth_masked_image,
-        landsat_ndvi_masked_image
+        modis_nbar_masked_image,
+        gch_masked_image,
+        s1_sar_masked_image,
+        gedi_agb_masked_image,
+        gedi_cth_masked_image,
+        landsat_ndvi_masked_image,
     ]
 
+
 @pytest.fixture
 def runner():
-    """ click runner for command line execution. """
+    """click runner for command line execution."""
     return CliRunner()
 
 
 @pytest.fixture
 def region_25ha_file() -> pathlib.Path:
-    """ Path to region_25ha geojson file. """
+    """Path to region_25ha geojson file."""
     return root_path.joinpath('tests/data/region_25ha.geojson')
 
 
 @pytest.fixture
 def region_100ha_file() -> pathlib.Path:
-    """ Path to region_100ha geojson file. """
+    """Path to region_100ha geojson file."""
     return root_path.joinpath('tests/data/region_100ha.geojson')
 
 
 @pytest.fixture
 def region_10000ha_file() -> pathlib.Path:
-    """ Path to region_10000ha geojson file. """
+    """Path to region_10000ha geojson file."""
     return root_path.joinpath('tests/data/region_10000ha.geojson')
 
 
@@ -343,4 +401,3 @@ def get_image_std(ee_image: ee.Image, region: Dict, std_scale: float):
         reducer='mean', geometry=region, crs=proj.crs(), scale=std_scale, bestEffort=True, maxPixels=1e6
     )
     return mean_std_image.get('TEST').getInfo()
-
diff --git a/tests/integration.py b/tests/integration.py
index aae678d..ffdff3d 100644
--- a/tests/integration.py
+++ b/tests/integration.py
@@ -13,31 +13,35 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
-import ee
-from pathlib import Path
-from httplib2 import Http
+
 import json
+from pathlib import Path
+
+import ee
 import numpy as np
+import pytest
 import rasterio as rio
+from click.testing import CliRunner
+from httplib2 import Http
+from rasterio.coords import BoundingBox
 from rasterio.crs import CRS
 from rasterio.features import bounds
 from rasterio.warp import transform_geom
-from rasterio.coords import BoundingBox
-from click.testing import CliRunner
-import pytest
+
 import geedim as gd
-from geedim import utils, cli
+from geedim import cli, utils
+from geedim.download import BaseImage
 
 
 @pytest.fixture(scope='session', autouse=True)
 def ee_init():
-    """ Override the ee_init fixture, so that we only initialise as geemap does, below. """
+    """Override the ee_init fixture, so that we only initialise as geemap does, below."""
     return
 
 
 def test_geemap_integration(tmp_path: Path):
-    """ Simulate the geemap download example. """
-    gd.Initialize(opt_url=None, http_transport=Http())    # a replica of geemap Initialize
+    """Simulate the geemap download example."""
+    gd.Initialize(opt_url=None, http_transport=Http())  # a replica of geemap Initialize
     ee_image = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA").first()
     gd_image = gd.download.BaseImage(ee_image)
     out_file = tmp_path.joinpath('landsat.tif')
@@ -47,27 +51,29 @@ def test_geemap_integration(tmp_path: Path):
 
 
 def test_geeml_integration(tmp_path: Path):
-    """ Test the geeml `user memory limit exceeded` example. """
+    """Test the geeml `user memory limit exceeded` example."""
     gd.Initialize()
     region = {
         'geodesic': False,
         'crs': {'type': 'name', 'properties': {'name': 'EPSG:4326'}},
         'type': 'Polygon',
-        'coordinates': [[
-            [6.030749828407996, 53.66867883985145],
-            [6.114742307473171, 53.66867883985145],
-            [6.114742307473171, 53.76381042843971],
-            [6.030749828407996, 53.76381042843971],
-            [6.030749828407996, 53.66867883985145]
-        ]]
+        'coordinates': [
+            [
+                [6.030749828407996, 53.66867883985145],
+                [6.114742307473171, 53.66867883985145],
+                [6.114742307473171, 53.76381042843971],
+                [6.030749828407996, 53.76381042843971],
+                [6.030749828407996, 53.66867883985145],
+            ]
+        ],
     }  # yapf: disable
 
     ee_image = (
-        ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED').
-        filterDate('2019-01-01', '2020-01-01').
-        filterBounds(region).
-        select(['B4', 'B3', 'B2', 'B8']).
-        reduce(ee.Reducer.percentile([35]))
+        ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
+        .filterDate('2019-01-01', '2020-01-01')
+        .filterBounds(region)
+        .select(['B4', 'B3', 'B2', 'B8'])
+        .reduce(ee.Reducer.percentile([35]))
     )  # yapf: disable
 
     gd_image = gd.download.BaseImage(ee_image)
@@ -82,26 +88,32 @@ def test_geeml_integration(tmp_path: Path):
 
     # test we can download the image with a max_tile_size of 16 MB
     gd_image.download(
-        out_file, crs='EPSG:4326', region=region, scale=10, dtype='float64', overwrite=True, max_tile_size=16,
+        out_file,
+        crs='EPSG:4326',
+        region=region,
+        scale=10,
+        dtype='float64',
+        overwrite=True,
+        max_tile_size=16,
     )
     assert out_file.exists()
     with rio.open(out_file, 'r') as ds:
         assert ds.count == 4
         assert ds.dtypes[0] == 'float64'
-        assert np.isnan(ds.nodata)
+        assert np.isinf(ds.nodata)
         region_cnrs = np.array(region['coordinates'][0])
         region_bounds = rio.coords.BoundingBox(*region_cnrs.min(axis=0), *region_cnrs.max(axis=0))
         # sometimes the top/bottom bounds of the dataset are swapped, so extract and compare UL and BR corners
         print(f'region_bounds: {region_bounds}')
         print(f'ds.bounds: {ds.bounds}')
-        ds_ul = np.array([min(ds.bounds.left,ds.bounds.right), min(ds.bounds.top,ds.bounds.bottom)])
-        ds_lr = np.array([max(ds.bounds.left,ds.bounds.right), max(ds.bounds.top,ds.bounds.bottom)])
+        ds_ul = np.array([min(ds.bounds.left, ds.bounds.right), min(ds.bounds.top, ds.bounds.bottom)])
+        ds_lr = np.array([max(ds.bounds.left, ds.bounds.right), max(ds.bounds.top, ds.bounds.bottom)])
         assert region_cnrs.min(axis=0) == pytest.approx(ds_ul, abs=1e-3)
         assert region_cnrs.max(axis=0) == pytest.approx(ds_lr, abs=1e-3)
 
 
 def test_cli_asset_export(l8_image_id, region_25ha_file: Path, runner: CliRunner, tmp_path: Path):
-    """  Export a test image to an asset using the CLI. """
+    """Export a test image to an asset using the CLI."""
     # create a randomly named folder to allow parallel tests without overwriting the same asset
     gd.Initialize()
     folder = f'geedim/int_test_asset_export_{np.random.randint(1 << 31)}'
@@ -118,7 +130,7 @@ def test_cli_asset_export(l8_image_id, region_25ha_file: Path, runner: CliRunner
             f'--dtype uint16 --mask --resampling bilinear --wait --type asset'
         )
         result = runner.invoke(cli.cli, cli_str.split())
-        assert (result.exit_code == 0)
+        assert result.exit_code == 0
         assert ee.data.getAsset(test_asset_id) is not None
 
         # download the asset image
@@ -139,13 +151,43 @@ def test_cli_asset_export(l8_image_id, region_25ha_file: Path, runner: CliRunner
     with open(region_25ha_file) as f:
         region = json.load(f)
     with rio.open(download_filename, 'r') as im:
-        im : rio.DatasetReader
+        im: rio.DatasetReader
         exp_region = transform_geom('EPSG:4326', im.crs, region)
         exp_bounds = BoundingBox(*bounds(exp_region))
         assert im.crs == CRS.from_string(crs)
         assert im.transform[0] == scale
         assert im.count > 1
         assert (
-            (im.bounds[0] <= exp_bounds[0]) and (im.bounds[1] <= exp_bounds[1]) and
-            (im.bounds[2] >= exp_bounds[2]) and (im.bounds[3] >= exp_bounds[3])
+            (im.bounds[0] <= exp_bounds[0])
+            and (im.bounds[1] <= exp_bounds[1])
+            and (im.bounds[2] >= exp_bounds[2])
+            and (im.bounds[3] >= exp_bounds[3])
         )
+
+
+@pytest.mark.parametrize('dtype', ['float32', 'float64', 'uint8', 'int8', 'uint16', 'int16', 'uint32', 'int32'])
+def test_ee_geotiff_nodata(dtype: str, l9_image_id: str):
+    """Test the nodata value of the Earth Engine GeoTIFF returned by ``ee.data.computePixels()`` or
+    ``ee.Image.getDownloadUrl()`` equals the geedim expected value (see
+    https://issuetracker.google.com/issues/350528377 for context).
+    """
+    # use geedim to prepare an image for downloading as dtype
+    gd.Initialize()
+    masked_image = gd.MaskedImage.from_id(l9_image_id)
+    shape = (10, 10)
+    exp_image, profile = masked_image._prepare_for_download(shape=shape, dtype=dtype)
+
+    # download a small tile with ee.data.computePixels
+    request = {
+        'expression': exp_image.ee_image,
+        'bandIds': ['SR_B3'],
+        'grid': {'dimensions': {'width': shape[1], 'height': shape[0]}},
+        'fileFormat': 'GEO_TIFF',
+    }
+    im_bytes = ee.data.computePixels(request)
+
+    # test nodata with rasterio
+    with rio.MemoryFile(im_bytes) as mf, mf.open() as ds:
+        assert ds.nodata == profile['nodata']
+        # test the EE dtype is not lower precision compared to expected dtype
+        assert np.promote_types(profile['dtype'], ds.dtypes[0]) == ds.dtypes[0]
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 66fa3c1..6a14773 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -13,62 +13,72 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
+
 import json
 import pathlib
 from datetime import datetime
 from glob import glob
-from typing import List, Dict, Tuple
+from typing import Dict, List, Tuple
 
 import ee
 import numpy as np
 import pytest
 import rasterio as rio
 from click.testing import CliRunner
-from geedim.cli import cli
-from geedim.utils import root_path, asset_id
 from rasterio.coords import BoundingBox
 from rasterio.crs import CRS
 from rasterio.features import bounds
-from rasterio.warp import transform_geom
 from rasterio.transform import Affine
+from rasterio.warp import transform_geom
+
+from geedim.cli import cli
+from geedim.utils import asset_id
 
 
 @pytest.fixture()
 def l4_5_image_id_list(l4_image_id, l5_image_id) -> List[str]:
-    """ A list of landsat 4 & 5 image ID's. """
+    """A list of landsat 4 & 5 image ID's."""
     return [l4_image_id, l5_image_id]
 
 
 @pytest.fixture()
 def l8_9_image_id_list(l8_image_id, l9_image_id) -> List[str]:
-    """ A list of landsat 8 & 9 image ID's. """
+    """A list of landsat 8 & 9 image ID's."""
     return [l8_image_id, l9_image_id]
 
 
 @pytest.fixture()
 def s2_sr_image_id_list() -> List[str]:
-    """ A list of Sentinel-2 SR image IDs. """
+    """A list of Sentinel-2 SR image IDs."""
     return [
         'COPERNICUS/S2_SR/20211004T080801_20211004T083709_T34HEJ',
         'COPERNICUS/S2_SR/20211123T081241_20211123T083704_T34HEJ',
-        'COPERNICUS/S2_SR/20220107T081229_20220107T083059_T34HEJ'
+        'COPERNICUS/S2_SR/20220107T081229_20220107T083059_T34HEJ',
     ]
 
 
 @pytest.fixture()
 def gedi_image_id_list() -> List[str]:
-    """ A list of GEDI canopy top height ID's. """
+    """A list of GEDI canopy top height ID's."""
     return [
-        'LARSE/GEDI/GEDI02_A_002_MONTHLY/202008_018E_036S', 'LARSE/GEDI/GEDI02_A_002_MONTHLY/202009_018E_036S',
-        'LARSE/GEDI/GEDI02_A_002_MONTHLY/202005_018E_036S'
+        'LARSE/GEDI/GEDI02_A_002_MONTHLY/202008_018E_036S',
+        'LARSE/GEDI/GEDI02_A_002_MONTHLY/202009_018E_036S',
+        'LARSE/GEDI/GEDI02_A_002_MONTHLY/202005_018E_036S',
     ]
 
 
 def _test_downloaded_file(
-    filename: pathlib.Path, region: Dict = None, crs: str = None, scale: float = None, dtype: str = None,
-    bands: List[str] = None, scale_offset: bool = None, transform: Affine = None, shape: Tuple[int, int] = None
+    filename: pathlib.Path,
+    region: Dict = None,
+    crs: str = None,
+    scale: float = None,
+    dtype: str = None,
+    bands: List[str] = None,
+    scale_offset: bool = None,
+    transform: Affine = None,
+    shape: Tuple[int, int] = None,
 ):
-    """ Helper function to test image file format against given parameters. """
+    """Helper function to test image file format against given parameters."""
     with rio.open(filename, 'r') as ds:
         ds: rio.DatasetReader = ds
         assert ds.nodata is not None
@@ -79,8 +89,10 @@ def _test_downloaded_file(
             exp_region = transform_geom('EPSG:4326', ds.crs, region)
             exp_bounds = BoundingBox(*bounds(exp_region))
             assert (
-                (ds.bounds[0] <= exp_bounds[0]) and (ds.bounds[1] <= exp_bounds[1]) and
-                (ds.bounds[2] >= exp_bounds[2]) and (ds.bounds[3] >= exp_bounds[3])
+                (ds.bounds[0] <= exp_bounds[0])
+                and (ds.bounds[1] <= exp_bounds[1])
+                and (ds.bounds[2] >= exp_bounds[2])
+                and (ds.bounds[3] >= exp_bounds[3])
             )
         if crs:
             assert CRS(ds.crs) == CRS.from_string(crs)
@@ -90,7 +102,8 @@ def _test_downloaded_file(
             assert ds.dtypes[0] == dtype
         if scale_offset:
             refl_bands = [
-                i for i in range(1, ds.count + 1)
+                i
+                for i in range(1, ds.count + 1)
                 if ('center_wavelength' in ds.tags(i)) and (float(ds.tags(i)['center_wavelength']) < 1)
             ]
             array = ds.read(refl_bands, masked=True)
@@ -102,11 +115,12 @@ def _test_downloaded_file(
             assert ds.shape == tuple(shape)
         if bands:
             assert set(bands) == set(ds.descriptions)
-            assert set(bands) ==  set([ds.tags(bi)['name'] for bi in range(1, ds.count + 1)])
+            assert set(bands) == set([ds.tags(bi)['name'] for bi in range(1, ds.count + 1)])
 
 
 @pytest.mark.parametrize(
-    'name, start_date, end_date, region, fill_portion, cloudless_portion, is_csmask', [
+    'name, start_date, end_date, region, fill_portion, cloudless_portion, is_csmask',
+    [
         ('LANDSAT/LC09/C02/T1_L2', '2022-01-01', '2022-02-01', 'region_100ha_file', 10, 50, True),
         ('LANDSAT/LE07/C02/T1_L2', '2022-01-01', '2022-02-01', 'region_100ha_file', 0, 0, True),
         ('LANDSAT/LT05/C02/T1_L2', '2005-01-01', '2006-02-01', 'region_100ha_file', 40, 50, True),
@@ -114,12 +128,20 @@ def _test_downloaded_file(
         ('COPERNICUS/S2', '2022-01-01', '2022-01-15', 'region_100ha_file', 50, 40, True),
         ('COPERNICUS/S2_SR_HARMONIZED', '2022-01-01', '2022-01-15', 'region_100ha_file', 0, 50, True),
         ('COPERNICUS/S2_HARMONIZED', '2022-01-01', '2022-01-15', 'region_100ha_file', 50, 40, True),
-        ('LARSE/GEDI/GEDI02_A_002_MONTHLY', '2021-11-01', '2022-01-01', 'region_100ha_file', 1, 0, False)
-    ]
+        ('LARSE/GEDI/GEDI02_A_002_MONTHLY', '2021-11-01', '2022-01-01', 'region_100ha_file', 1, 0, False),
+    ],
 )
 def test_search(
-    name, start_date: str, end_date: str, region: str, fill_portion: float, cloudless_portion: float, is_csmask: bool,
-    tmp_path: pathlib.Path, runner: CliRunner, request: pytest.FixtureRequest
+    name,
+    start_date: str,
+    end_date: str,
+    region: str,
+    fill_portion: float,
+    cloudless_portion: float,
+    is_csmask: bool,
+    tmp_path: pathlib.Path,
+    runner: CliRunner,
+    request: pytest.FixtureRequest,
 ):
     """
     Test search command gives valid results for different cloud/shadow maskable, and generic collections.
@@ -131,8 +153,8 @@ def test_search(
         f'-cp {cloudless_portion} -op {results_file}'
     )
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
-    assert (results_file.exists())
+    assert result.exit_code == 0
+    assert results_file.exists()
     with open(results_file, 'r') as f:
         properties = json.load(f)
 
@@ -155,7 +177,7 @@ def test_search(
 
 
 def test_config_search_s2(region_10000ha_file: pathlib.Path, runner: CliRunner, tmp_path: pathlib.Path):
-    """ Test `config` sub-command chained with `search` of Sentinel-2 affects CLOUDLESS_PORTION as expected. """
+    """Test `config` sub-command chained with `search` of Sentinel-2 affects CLOUDLESS_PORTION as expected."""
     results_file = tmp_path.joinpath('search_results.json')
     name = 'COPERNICUS/S2_SR'
     cl_portion_list = []
@@ -165,8 +187,8 @@ def test_config_search_s2(region_10000ha_file: pathlib.Path, runner: CliRunner,
             f'{results_file}'
         )
         result = runner.invoke(cli, cli_str.split())
-        assert (result.exit_code == 0)
-        assert (results_file.exists())
+        assert result.exit_code == 0
+        assert results_file.exists()
         with open(results_file, 'r') as f:
             properties = json.load(f)
         cl_portion_list.append(np.array([prop_dict['CLOUDLESS_PORTION'] for prop_dict in properties.values()]))
@@ -176,7 +198,7 @@ def test_config_search_s2(region_10000ha_file: pathlib.Path, runner: CliRunner,
 
 
 def test_config_search_l9(region_10000ha_file: pathlib.Path, runner: CliRunner, tmp_path: pathlib.Path):
-    """ Test `config` sub-command chained with `search` of Landsat-9 affects CLOUDLESS_PORTION as expected. """
+    """Test `config` sub-command chained with `search` of Landsat-9 affects CLOUDLESS_PORTION as expected."""
     results_file = tmp_path.joinpath('search_results.json')
     name = 'LANDSAT/LC09/C02/T1_L2'
     cl_portion_list = []
@@ -186,8 +208,8 @@ def test_config_search_l9(region_10000ha_file: pathlib.Path, runner: CliRunner,
             f' {results_file}'
         )
         result = runner.invoke(cli, cli_str.split())
-        assert (result.exit_code == 0)
-        assert (results_file.exists())
+        assert result.exit_code == 0
+        assert results_file.exists()
         with open(results_file, 'r') as f:
             properties = json.load(f)
         cl_portion_list.append(np.array([prop_dict['CLOUDLESS_PORTION'] for prop_dict in properties.values()]))
@@ -197,7 +219,7 @@ def test_config_search_l9(region_10000ha_file: pathlib.Path, runner: CliRunner,
 
 
 def test_region_bbox_search(region_100ha_file: pathlib.Path, runner: CliRunner, tmp_path: pathlib.Path):
-    """ Test --bbox gives same search results as --region <geojson file>. """
+    """Test --bbox gives same search results as --region <geojson file>."""
 
     results_file = tmp_path.joinpath('search_results.json')
     with open(region_100ha_file, 'r') as f:
@@ -206,14 +228,14 @@ def test_region_bbox_search(region_100ha_file: pathlib.Path, runner: CliRunner,
     bbox_str = ' '.join([str(b) for b in bbox])
     cli_strs = [
         f'search -c LANDSAT/LC09/C02/T1_L2 -s 2022-01-01 -e 2022-02-01 -r {region_100ha_file} -op {results_file}',
-        f'search -c LANDSAT/LC09/C02/T1_L2 -s 2022-01-01 -e 2022-02-01 -b {bbox_str} -op {results_file}'
+        f'search -c LANDSAT/LC09/C02/T1_L2 -s 2022-01-01 -e 2022-02-01 -b {bbox_str} -op {results_file}',
     ]
 
     props_list = []
     for cli_str in cli_strs:
         result = runner.invoke(cli, cli_str.split())
-        assert (result.exit_code == 0)
-        assert (results_file.exists())
+        assert result.exit_code == 0
+        assert results_file.exists()
         with open(results_file, 'r') as f:
             properties = json.load(f)
         props_list.append(properties)
@@ -222,19 +244,19 @@ def test_region_bbox_search(region_100ha_file: pathlib.Path, runner: CliRunner,
 
 
 def test_raster_region_search(const_image_25ha_file, region_25ha_file, runner: CliRunner, tmp_path: pathlib.Path):
-    """ Test --region works with a raster file. """
+    """Test --region works with a raster file."""
 
     results_file = tmp_path.joinpath('search_results.json')
     cli_strs = [
         f'search -c LANDSAT/LC09/C02/T1_L2 -s 2022-01-01 -e 2022-02-01 -r {region_25ha_file} -op {results_file}',
-        f'search -c LANDSAT/LC09/C02/T1_L2 -s 2022-01-01 -e 2022-02-01 -r {const_image_25ha_file} -op {results_file}'
+        f'search -c LANDSAT/LC09/C02/T1_L2 -s 2022-01-01 -e 2022-02-01 -r {const_image_25ha_file} -op {results_file}',
     ]
 
     props_list = []
     for cli_str in cli_strs:
         result = runner.invoke(cli, cli_str.split())
-        assert (result.exit_code == 0)
-        assert (results_file.exists())
+        assert result.exit_code == 0
+        assert results_file.exists()
         with open(results_file, 'r') as f:
             properties = json.load(f)
         props_list.append(properties)
@@ -243,17 +265,15 @@ def test_raster_region_search(const_image_25ha_file, region_25ha_file, runner: C
 
 
 def test_search_add_props_l9(region_25ha_file: pathlib.Path, runner: CliRunner, tmp_path: pathlib.Path):
-    """ Test --add-property generates results with the additional property keys. """
+    """Test --add-property generates results with the additional property keys."""
     results_file = tmp_path.joinpath('search_results.json')
     name = 'LANDSAT/LC09/C02/T1_L2'
     add_props = ['CLOUD_COVER', 'GEOMETRIC_RMSE_VERIFY']
     add_props_str = ''.join([f' -ap {add_prop} ' for add_prop in add_props])
-    cli_str = (
-        f'search -c {name} -s 2022-01-01 -e 2022-02-01 -r {region_25ha_file} {add_props_str} -op {results_file}'
-    )
+    cli_str = f'search -c {name} -s 2022-01-01 -e 2022-02-01 -r {region_25ha_file} {add_props_str} -op {results_file}'
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
-    assert (results_file.exists())
+    assert result.exit_code == 0
+    assert results_file.exists()
     with open(results_file, 'r') as f:
         properties = json.load(f)
     prop_keys = list(properties.values())[0].keys()
@@ -262,7 +282,7 @@ def test_search_add_props_l9(region_25ha_file: pathlib.Path, runner: CliRunner,
 
 
 def test_search_custom_filter_l9(region_25ha_file: pathlib.Path, runner: CliRunner, tmp_path: pathlib.Path):
-    """ Test --custom-filter filters the search results as specified. """
+    """Test --custom-filter filters the search results as specified."""
     results_file = tmp_path.joinpath('search_results.json')
     name = 'LANDSAT/LC09/C02/T1_L2'
     cc_thresh = 30
@@ -272,8 +292,8 @@ def test_search_custom_filter_l9(region_25ha_file: pathlib.Path, runner: CliRunn
         f'-cf {add_prop}>{cc_thresh} -op {results_file}'
     )
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
-    assert (results_file.exists())
+    assert result.exit_code == 0
+    assert results_file.exists()
     with open(results_file, 'r') as f:
         properties = json.load(f)
     prop_keys = list(properties.values())[0].keys()
@@ -281,20 +301,16 @@ def test_search_custom_filter_l9(region_25ha_file: pathlib.Path, runner: CliRunn
     assert all([prop[add_prop] > cc_thresh for prop in properties.values()])
 
 
-def test_search_cloudless_portion_no_value(
-    region_25ha_file: pathlib.Path, runner: CliRunner, tmp_path: pathlib.Path
-):
-    """ Test `search --cloudless-portion` gives the same results as `search --cloudless-portion 0`. """
+def test_search_cloudless_portion_no_value(region_25ha_file: pathlib.Path, runner: CliRunner, tmp_path: pathlib.Path):
+    """Test `search --cloudless-portion` gives the same results as `search --cloudless-portion 0`."""
     results_file = tmp_path.joinpath(f'search_results.json')
     name = 'LANDSAT/LC09/C02/T1_L2'
     clp_list = []
     for post_fix, cp_spec in zip(['no_val', 'zero_val'], ['-cp', '-cp 0']):
-        cli_str = (
-            f'search -c {name} -s 2022-01-01 -e 2022-04-01 -r {region_25ha_file} {cp_spec} -op {results_file}'
-        )
+        cli_str = f'search -c {name} -s 2022-01-01 -e 2022-04-01 -r {region_25ha_file} {cp_spec} -op {results_file}'
         result = runner.invoke(cli, cli_str.split())
-        assert (result.exit_code == 0)
-        assert (results_file.exists())
+        assert result.exit_code == 0
+        assert results_file.exists()
         with open(results_file, 'r') as f:
             props = json.load(f)
         prop_keys = list(props.values())[0].keys()
@@ -306,25 +322,26 @@ def test_search_cloudless_portion_no_value(
 
 
 @pytest.mark.parametrize(
-    'image_id, region_file', [
+    'image_id, region_file',
+    [
         ('l8_image_id', 'region_25ha_file'),
         ('s2_sr_hm_image_id', 'region_25ha_file'),
         ('gedi_cth_image_id', 'region_25ha_file'),
         ('modis_nbar_image_id', 'region_25ha_file'),
-    ]
+    ],
 )  # yapf: disable
 def test_download_region_defaults(
     image_id: str, region_file: pathlib.Path, tmp_path: pathlib.Path, runner: CliRunner, request
 ):
-    """ Test image download with default crs, scale, dtype etc.  """
+    """Test image download with default crs, scale, dtype etc."""
     image_id = request.getfixturevalue(image_id)
     region_file = request.getfixturevalue(region_file)
     out_file = tmp_path.joinpath(image_id.replace('/', '-') + '.tif')
 
     cli_str = f'download -i {image_id} -r {region_file} -dd {tmp_path}'
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
-    assert (out_file.exists())
+    assert result.exit_code == 0
+    assert out_file.exists()
 
     # test downloaded file readability and format
     with open(region_file) as f:
@@ -333,16 +350,17 @@ def test_download_region_defaults(
 
 
 @pytest.mark.parametrize(
-    'image_id, region_file', [
+    'image_id, region_file',
+    [
         ('l8_image_id', 'region_25ha_file'),
         ('s2_sr_hm_image_id', 'region_25ha_file'),
         ('gedi_cth_image_id', 'region_25ha_file'),
-    ]
+    ],
 )  # yapf: disable
 def test_download_crs_transform(
     image_id: str, region_file: pathlib.Path, tmp_path: pathlib.Path, runner: CliRunner, request
 ):
-    """ Test image download with crs, crs_transform, & shape specified. """
+    """Test image download with crs, crs_transform, & shape specified."""
     image_id = request.getfixturevalue(image_id)
     region_file = request.getfixturevalue(region_file)
     out_file = tmp_path.joinpath(image_id.replace('/', '-') + '.tif')
@@ -360,8 +378,8 @@ def test_download_crs_transform(
     # run the download
     cli_str = f'download -i {image_id} -c {crs} -ct {crs_transform_str} -sh {shape_str} -dd {tmp_path}'
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
-    assert (out_file.exists())
+    assert result.exit_code == 0
+    assert out_file.exists()
 
     # test downloaded file readability and format
     _test_downloaded_file(out_file, crs=crs, region=region, transform=crs_transform, shape=shape)
@@ -370,21 +388,22 @@ def test_download_crs_transform(
 def test_download_like(
     l8_image_id: str, s2_sr_image_id: str, region_25ha_file: pathlib.Path, tmp_path: pathlib.Path, runner: CliRunner
 ):
-    """ Test image download using --like. """
+    """Test image download using --like."""
     l8_file = tmp_path.joinpath(l8_image_id.replace('/', '-') + '.tif')
     s2_file = tmp_path.joinpath(s2_sr_image_id.replace('/', '-') + '.tif')
 
     # download the landsat 8 image to be the template
+    # TODO: make a synthetic template to save time
     l8_cli_str = f'download -i {l8_image_id} -r {region_25ha_file} -dd {tmp_path}'
     result = runner.invoke(cli, l8_cli_str.split())
-    assert (result.exit_code == 0)
-    assert (l8_file.exists())
+    assert result.exit_code == 0
+    assert l8_file.exists()
 
     # download the sentinel 2 image like the landsat 8 image
     s2_cli_str = f'download -i {s2_sr_image_id} --like {l8_file} -dd {tmp_path}'
     result = runner.invoke(cli, s2_cli_str.split())
-    assert (result.exit_code == 0)
-    assert (s2_file.exists())
+    assert result.exit_code == 0
+    assert s2_file.exists()
 
     # test the landsat 8 image is 'like' the sentinel 2 image
     with rio.open(l8_file) as l8_im, rio.open(s2_file, 'r') as s2_im:
@@ -394,27 +413,57 @@ def test_download_like(
 
 
 @pytest.mark.parametrize(
-    'image_id, region_file, crs, scale, dtype, bands, mask, resampling, scale_offset, max_tile_size, max_tile_dim', [
+    'image_id, region_file, crs, scale, dtype, bands, mask, resampling, scale_offset, max_tile_size, max_tile_dim',
+    [
         ('l5_image_id', 'region_25ha_file', 'EPSG:3857', 30, 'uint16', None, False, 'near', False, 16, 10000),
         ('l9_image_id', 'region_25ha_file', 'EPSG:3857', 30, 'float32', None, False, 'near', True, 32, 10000),
         (
-            's2_toa_image_id', 'region_25ha_file', 'EPSG:3857', 10, 'float64', ['B5', 'B9'], True, 'bilinear', True, 32,
-            10000
+            's2_toa_image_id',
+            'region_25ha_file',
+            'EPSG:3857',
+            10,
+            'float64',
+            ['B5', 'B9'],
+            True,
+            'bilinear',
+            True,
+            32,
+            10000,
         ),
         ('modis_nbar_image_id', 'region_100ha_file', 'EPSG:3857', 500, 'int32', None, False, 'bicubic', False, 4, 100),
         (
-            'gedi_cth_image_id', 'region_25ha_file', 'EPSG:3857', 10, 'float32', ['rh99'], True, 'bilinear', False, 32,
-            10000
+            'gedi_cth_image_id',
+            'region_25ha_file',
+            'EPSG:3857',
+            10,
+            'float32',
+            ['rh99'],
+            True,
+            'bilinear',
+            False,
+            32,
+            10000,
         ),
         ('landsat_ndvi_image_id', 'region_25ha_file', 'EPSG:3857', 30, 'float64', None, True, 'near', False, 32, 10000),
-    ]
-) # yapf: disable
+    ],
+)  # yapf: disable
 def test_download_params(
-    image_id: str, region_file: str, crs: str, scale: float, dtype: str, bands: List[str], mask: bool, resampling: str,
-    scale_offset: bool, max_tile_size: float, max_tile_dim: int, tmp_path: pathlib.Path, runner: CliRunner,
-    request: pytest.FixtureRequest
+    image_id: str,
+    region_file: str,
+    crs: str,
+    scale: float,
+    dtype: str,
+    bands: List[str],
+    mask: bool,
+    resampling: str,
+    scale_offset: bool,
+    max_tile_size: float,
+    max_tile_dim: int,
+    tmp_path: pathlib.Path,
+    runner: CliRunner,
+    request: pytest.FixtureRequest,
 ):
-    """ Test image download, specifying all cli params except crs_transform and shape. """
+    """Test image download, specifying all cli params except crs_transform and shape."""
     image_id = request.getfixturevalue(image_id)
     region_file = request.getfixturevalue(region_file)
     out_file = tmp_path.joinpath(image_id.replace('/', '-') + '.tif')
@@ -427,8 +476,8 @@ def test_download_params(
     cli_str += ' --scale-offset' if scale_offset else ' --no-scale-offset'
     cli_str += ''.join([f' --band-name {bn}' for bn in bands]) if bands else ''
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
-    assert (out_file.exists())
+    assert result.exit_code == 0
+    assert out_file.exists()
 
     with open(region_file) as f:
         region = json.load(f)
@@ -441,10 +490,10 @@ def test_download_params(
 def test_max_tile_size_error(
     s2_sr_image_id: str, region_100ha_file: pathlib.Path, tmp_path: pathlib.Path, runner: CliRunner, request
 ):
-    """ Test image download with max_tile_size > EE limit raises an EE error.  """
+    """Test image download with max_tile_size > EE limit raises an EE error."""
     cli_str = f'download -i {s2_sr_image_id} -r {region_100ha_file} -dd {tmp_path} -mts 100'
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code != 0)
+    assert result.exit_code != 0
     assert isinstance(result.exception, ValueError)
     assert 'download size limit' in str(result.exception)
 
@@ -452,26 +501,26 @@ def test_max_tile_size_error(
 def test_max_tile_dim_error(
     s2_sr_image_id: str, region_100ha_file: pathlib.Path, tmp_path: pathlib.Path, runner: CliRunner, request
 ):
-    """ Test image download with max_tile_dim > EE limit raises an EE error.  """
+    """Test image download with max_tile_dim > EE limit raises an EE error."""
     cli_str = f'download -i {s2_sr_image_id} -r {region_100ha_file} -dd {tmp_path} -mtd 100000'
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code != 0)
+    assert result.exit_code != 0
     assert isinstance(result.exception, ValueError)
     assert 'download limit' in str(result.exception)
 
 
 def test_export_drive_params(l8_image_id: str, region_25ha_file: pathlib.Path, runner: CliRunner):
-    """ Test export to google drive starts ok, specifying all cli params"""
+    """Test export to google drive starts ok, specifying all cli params"""
     cli_str = (
         f'export -i {l8_image_id} -r {region_25ha_file} -df geedim/test --crs EPSG:3857 --scale 30 --dtype uint16 '
         f'--mask --resampling bilinear --no-wait --band-name SR_B4 --band-name SR_B3 --band-name SR_B2'
     )
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
+    assert result.exit_code == 0
 
 
 def test_export_asset_params(l8_image_id: str, region_25ha_file: pathlib.Path, runner: CliRunner):
-    """ Test export to asset starts ok, specifying all cli params"""
+    """Test export to asset starts ok, specifying all cli params"""
     # Note when e.g. github runs this test in parallel, it could run into problems trying to overwrite an existing
     # asset.  The overwrite error won't be raised with --no-wait though.  So this test serves at least to check the
     # CLI export options work, and won't fail if run in parallel, even if it runs into overwrite problems.
@@ -487,11 +536,11 @@ def test_export_asset_params(l8_image_id: str, region_25ha_file: pathlib.Path, r
         f'--mask --resampling bilinear --no-wait --type asset --band-name SR_B4 --band-name SR_B3 --band-name SR_B2'
     )
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
+    assert result.exit_code == 0
 
 
 def test_export_asset_no_folder_error(l8_image_id: str, region_25ha_file: pathlib.Path, runner: CliRunner):
-    """ Test export to asset raises an error when no folder is specified. """
+    """Test export to asset raises an error when no folder is specified."""
     cli_str = (
         f'export -i {l8_image_id} -r {region_25ha_file} --crs EPSG:3857 --scale 30 '
         f'--dtype uint16 --mask --resampling bilinear --no-wait --type asset'
@@ -503,16 +552,22 @@ def test_export_asset_no_folder_error(l8_image_id: str, region_25ha_file: pathli
 
 @pytest.mark.parametrize('image_list, scale', [('s2_sr_image_id_list', 10), ('l8_9_image_id_list', 30)])
 def test_composite_defaults(
-    image_list: str, scale: float, region_25ha_file: pathlib.Path, runner: CliRunner, tmp_path: pathlib.Path,
-    request: pytest.FixtureRequest
+    image_list: str,
+    scale: float,
+    region_25ha_file: pathlib.Path,
+    runner: CliRunner,
+    tmp_path: pathlib.Path,
+    request: pytest.FixtureRequest,
 ):
-    """ Test composite with default CLI parameters.  """
+    """Test composite with default CLI parameters."""
     image_list = request.getfixturevalue(image_list)
     image_ids_str = ' -i '.join(image_list)
-    cli_str = f'composite -i {image_ids_str} download --crs EPSG:3857 --scale {scale} -r {region_25ha_file} -dd' \
-              f' {tmp_path}'
+    cli_str = (
+        f'composite -i {image_ids_str} download --crs EPSG:3857 --scale {scale} -r {region_25ha_file} -dd'
+        f' {tmp_path}'
+    )
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
+    assert result.exit_code == 0
 
     # test downloaded file exists
     out_files = glob(str(tmp_path.joinpath(f'*COMP*.tif')))
@@ -525,18 +580,28 @@ def test_composite_defaults(
 
 
 @pytest.mark.parametrize(
-    'image_list, method, region_file, date, mask, resampling, download_scale', [
+    'image_list, method, region_file, date, mask, resampling, download_scale',
+    [
         ('s2_sr_image_id_list', 'mosaic', None, '2021-10-01', True, 'near', 10),
         ('l8_9_image_id_list', 'q-mosaic', 'region_25ha_file', None, True, 'bilinear', 30),
         ('l8_9_image_id_list', 'medoid', 'region_25ha_file', None, True, 'near', 30),
         ('gedi_image_id_list', 'medoid', None, None, True, 'bilinear', 25),
-    ]
+    ],
 )
 def test_composite_params(
-    image_list: str, method: str, region_file: str, date: str, mask: bool, resampling: str, download_scale: float,
-    region_25ha_file, runner: CliRunner, tmp_path: pathlib.Path, request: pytest.FixtureRequest
+    image_list: str,
+    method: str,
+    region_file: str,
+    date: str,
+    mask: bool,
+    resampling: str,
+    download_scale: float,
+    region_25ha_file,
+    runner: CliRunner,
+    tmp_path: pathlib.Path,
+    request: pytest.FixtureRequest,
 ):
-    """ Test composite with default CLI parameters. """
+    """Test composite with default CLI parameters."""
     image_list = request.getfixturevalue(image_list)
     region_file = request.getfixturevalue(region_file) if region_file else None
     image_ids_str = ' -i '.join(image_list)
@@ -547,7 +612,7 @@ def test_composite_params(
     cli_download_str = f'download -r {region_25ha_file} --crs EPSG:3857 --scale {download_scale} -dd {tmp_path}'
     cli_str = cli_comp_str + ' ' + cli_download_str
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
+    assert result.exit_code == 0
 
     # test downloaded file exists
     out_files = glob(str(tmp_path.joinpath(f'*COMP*.tif')))
@@ -560,14 +625,14 @@ def test_composite_params(
 
 
 def test_search_composite_download(region_25ha_file, runner: CliRunner, tmp_path: pathlib.Path):
-    """ Test chaining of `search`, `composite` and `download`. """
+    """Test chaining of `search`, `composite` and `download`."""
 
     cli_search_str = f'search -c COPERNICUS/S1_GRD -s 2022-01-01 -e 2022-02-01 -r {region_25ha_file}'
     cli_comp_str = f'composite --mask'
     cli_download_str = f'download --crs EPSG:3857 --scale 10 -dd {tmp_path}'
     cli_str = cli_search_str + ' ' + cli_comp_str + ' ' + cli_download_str
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
+    assert result.exit_code == 0
 
     # test downloaded file exists
     out_files = glob(str(tmp_path.joinpath(f'*COMP*.tif')))
@@ -591,7 +656,7 @@ def test_search_composite_x2_download(region_25ha_file, runner: CliRunner, tmp_p
     cli_download_str = f'download --crs EPSG:3857 --scale 30 -dd {tmp_path}'
     cli_str = cli_search_str + ' ' + cli_comp1_str + ' ' + cli_comp2_str + ' ' + cli_download_str
     result = runner.invoke(cli, cli_str.split())
-    assert (result.exit_code == 0)
+    assert result.exit_code == 0
 
     # test downloaded file exists
     out_files = glob(str(tmp_path.joinpath(f'*COMP*.tif')))
diff --git a/tests/test_collection.py b/tests/test_collection.py
index 7b04ee0..2071bde 100644
--- a/tests/test_collection.py
+++ b/tests/test_collection.py
@@ -13,55 +13,57 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
+
 from datetime import datetime
-from typing import List, Union, Dict
+from typing import Dict, List, Union
 
 import ee
 import numpy as np
 import pytest
+
 from geedim import schema
 from geedim.collection import MaskedCollection
 from geedim.enums import CompositeMethod, ResamplingMethod
-from geedim.errors import UnfilteredError, InputImageError
+from geedim.errors import InputImageError, UnfilteredError
 from geedim.mask import MaskedImage
-from geedim.utils import split_id, get_projection
-
+from geedim.utils import get_projection, split_id
 from .conftest import get_image_std
 
 
 @pytest.fixture()
 def l4_5_image_list(l4_image_id, l5_masked_image) -> List[Union[str, MaskedImage]]:
-    """ A list of landsat 4 & 5 image ID / MaskedImage's """
+    """A list of landsat 4 & 5 image ID / MaskedImage's"""
     return [l4_image_id, l5_masked_image]
 
 
 @pytest.fixture()
 def l8_9_image_list(l8_image_id, l9_masked_image) -> List[Union[str, MaskedImage]]:
-    """ A list of landsat 8 & 9 image IDs/ MaskedImage's """
+    """A list of landsat 8 & 9 image IDs/ MaskedImage's"""
     return [l8_image_id, l9_masked_image]
 
 
 @pytest.fixture()
 def s2_sr_image_list() -> List[Union[str, MaskedImage]]:
-    """ A list of Sentinel-2 SR image IDs/ MaskedImage's """
+    """A list of Sentinel-2 SR image IDs/ MaskedImage's"""
     return [
         'COPERNICUS/S2_SR/20211004T080801_20211004T083709_T34HEJ',
         'COPERNICUS/S2_SR/20211123T081241_20211123T083704_T34HEJ',
-        MaskedImage.from_id('COPERNICUS/S2_SR/20220107T081229_20220107T083059_T34HEJ')
+        MaskedImage.from_id('COPERNICUS/S2_SR/20220107T081229_20220107T083059_T34HEJ'),
     ]
 
 
 @pytest.fixture()
 def gedi_image_list() -> List[Union[str, MaskedImage]]:
-    """ A list of GEDI canopy top height IDs/ MaskedImage's """
+    """A list of GEDI canopy top height IDs/ MaskedImage's"""
     return [
-        'LARSE/GEDI/GEDI02_A_002_MONTHLY/202008_018E_036S', 'LARSE/GEDI/GEDI02_A_002_MONTHLY/202009_018E_036S',
-        MaskedImage.from_id('LARSE/GEDI/GEDI02_A_002_MONTHLY/202005_018E_036S')
+        'LARSE/GEDI/GEDI02_A_002_MONTHLY/202008_018E_036S',
+        'LARSE/GEDI/GEDI02_A_002_MONTHLY/202009_018E_036S',
+        MaskedImage.from_id('LARSE/GEDI/GEDI02_A_002_MONTHLY/202005_018E_036S'),
     ]
 
 
 def test_split_id():
-    """ Test split_id(). """
+    """Test split_id()."""
     coll_name, im_id = split_id('A/B/C')
     assert coll_name == 'A/B'
     assert im_id == 'C'
@@ -72,7 +74,7 @@ def test_split_id():
 
 @pytest.mark.parametrize('name', ['l9_image_id', 'gch_image_id'])
 def test_from_name(name: str, request):
-    """ Test MaskedCollection.from_name() for non Sentinel-2 collections. """
+    """Test MaskedCollection.from_name() for non Sentinel-2 collections."""
     name = request.getfixturevalue(name)
     name, _ = split_id(name)
     gd_collection = MaskedCollection.from_name(name)
@@ -105,7 +107,7 @@ def test_from_name_s2(name: str, request):
 
 
 def test_unfiltered_error(s2_sr_image_id):
-    """ Test UnfilteredError is raised when calling `properties` or `composite` on an unfiltered collection. """
+    """Test UnfilteredError is raised when calling `properties` or `composite` on an unfiltered collection."""
     gd_collection = MaskedCollection.from_name(split_id(s2_sr_image_id)[0])
     with pytest.raises(UnfilteredError):
         _ = gd_collection.properties
@@ -113,8 +115,8 @@ def test_unfiltered_error(s2_sr_image_id):
         _ = gd_collection.composite()
 
 
-def test_from_list_errors(landsat_image_ids, s2_image_ids, user_masked_image):
-    """ Test MaskedCollection.from_list() various error cases. """
+def test_from_list_errors(landsat_image_ids, user_masked_image):
+    """Test MaskedCollection.from_list() various error cases."""
     with pytest.raises(InputImageError):
         # test an error is raised when an image has no 'id'/'system:time_start' property
         MaskedCollection.from_list([landsat_image_ids[0], user_masked_image])
@@ -166,7 +168,7 @@ def test_from_list(image_list: str, request):
 
 @pytest.mark.parametrize('image_list', ['s2_sr_image_list', 'gedi_image_list'])
 def test_from_list_order(image_list: str, request):
-    """ Test MaskedCollection.from_list() maintains the order of the provided image list. """
+    """Test MaskedCollection.from_list() maintains the order of the provided image list."""
     image_list: List = request.getfixturevalue(image_list)[::-1]
     image_ids = [im_obj if isinstance(im_obj, str) else im_obj.id for im_obj in image_list]
     gd_collection = MaskedCollection.from_list(image_list)
@@ -174,19 +176,21 @@ def test_from_list_order(image_list: str, request):
 
 
 def test_from_list_ee_image(gedi_image_list: List):
-    """ Test MaskedCollection.from_list() with an ee.Image in the list. """
+    """Test MaskedCollection.from_list() with an ee.Image in the list."""
     image_ids = [im_obj if isinstance(im_obj, str) else im_obj.id for im_obj in gedi_image_list]
     image_list = gedi_image_list
     image_list[1] = ee.Image(image_list[1])
     gd_collection = MaskedCollection.from_list(image_list)
     assert list(gd_collection.properties.keys()) == image_ids
 
+
 @pytest.mark.parametrize(
-    'image_list, add_props', [
+    'image_list, add_props',
+    [
         ('s2_sr_image_list', ['AOT_RETRIEVAL_ACCURACY', 'CLOUDY_PIXEL_PERCENTAGE']),
-        ('l8_9_image_list', ['CLOUD_COVER', 'GEOMETRIC_RMSE_VERIFY'])
-    ]
-)  # yapf: disable
+        ('l8_9_image_list', ['CLOUD_COVER', 'GEOMETRIC_RMSE_VERIFY']),
+    ],
+)
 def test_from_list_add_props(image_list: str, add_props: List, request: pytest.FixtureRequest):
     """
     Test MaskedCollection.from_list(add_props=...) contains the add_props in properties and schema.
@@ -204,14 +208,15 @@ def test_from_list_add_props(image_list: str, add_props: List, request: pytest.F
 
 
 @pytest.mark.parametrize(
-    'name, start_date, end_date, region, fill_portion, cloudless_portion, is_csmask', [
+    'name, start_date, end_date, region, fill_portion, cloudless_portion, is_csmask',
+    [
         ('LANDSAT/LC09/C02/T1_L2', '2022-01-01', '2022-02-01', 'region_100ha', 0, 50, True),
         ('LANDSAT/LE07/C02/T1_L2', '2022-01-01', '2022-02-01', 'region_100ha', 0, 0, True),
         ('LANDSAT/LT05/C02/T1_L2', '2005-01-01', '2006-02-01', 'region_100ha', 40, 50, True),
         ('COPERNICUS/S2_SR', '2022-01-01', '2022-01-15', 'region_100ha', 0, 50, True),
         ('COPERNICUS/S2_HARMONIZED', '2022-01-01', '2022-01-15', 'region_100ha', 50, 40, True),
-        ('LARSE/GEDI/GEDI02_A_002_MONTHLY', '2021-08-01', '2021-09-01', 'region_100ha', .1, 0, False)
-    ]
+        ('LARSE/GEDI/GEDI02_A_002_MONTHLY', '2021-08-01', '2021-09-01', 'region_100ha', 0.1, 0, False),
+    ],
 )
 def test_search(
     name, start_date: str, end_date: str, region: str, fill_portion: float, cloudless_portion: float, is_csmask, request
@@ -247,7 +252,7 @@ def test_search(
 
 
 def test_empty_search(region_100ha):
-    """ Test MaskedCollection.search() for empty search results. """
+    """Test MaskedCollection.search() for empty search results."""
     gd_collection = MaskedCollection.from_name('LANDSAT/LC09/C02/T1_L2')
     searched_collection = gd_collection.search('2022-01-01', '2022-01-02', region_100ha, cloudless_portion=100)
     assert searched_collection.properties is not None
@@ -256,7 +261,7 @@ def test_empty_search(region_100ha):
 
 
 def test_search_no_end_date(region_100ha):
-    """ Test MaskedCollection.search() with ``end_date=None`` searches for a day from ``start_date`. """
+    """Test MaskedCollection.search() with ``end_date=None`` searches for a day from ``start_date``."""
     gd_collection = MaskedCollection.from_name('LANDSAT/LC09/C02/T1_L2')
     searched_collection = gd_collection.search('2022-01-03', None, region_100ha)
 
@@ -271,7 +276,7 @@ def test_search_no_end_date(region_100ha):
 
 
 def test_search_date_error(region_100ha):
-    """ Test MaskedCollection.search() raises an error when end date is on or before start date. """
+    """Test MaskedCollection.search() raises an error when end date is on or before start date."""
     gd_collection = MaskedCollection.from_name('LANDSAT/LC09/C02/T1_L2')
     with pytest.raises(ValueError):
         _ = gd_collection.search('2022-01-02', '2022-01-01', region_100ha)
@@ -328,11 +333,12 @@ def test_search_add_props(region_25ha):
 
 
 @pytest.mark.parametrize(
-    'name, start_date, end_date, region', [
+    'name, start_date, end_date, region',
+    [
         ('LANDSAT/LC09/C02/T1_L2', '2022-01-01', '2022-02-01', 'region_100ha'),
         ('COPERNICUS/S2_HARMONIZED', '2022-01-01', '2022-01-15', 'region_100ha'),
-        ('LARSE/GEDI/GEDI02_A_002_MONTHLY', '2021-08-01', '2021-09-01', 'region_100ha')
-    ]
+        ('LARSE/GEDI/GEDI02_A_002_MONTHLY', '2021-08-01', '2021-09-01', 'region_100ha'),
+    ],
 )
 def test_search_no_fill_or_cloudless_portion(
     name: str, start_date: str, end_date: str, region: str, request: pytest.FixtureRequest
@@ -363,14 +369,15 @@ def test_search_no_fill_or_cloudless_portion(
 
 
 @pytest.mark.parametrize(
-    'image_list, method, region, date', [
+    'image_list, method, region, date',
+    [
         ('s2_sr_image_list', CompositeMethod.q_mosaic, 'region_10000ha', None),
         ('s2_sr_image_list', CompositeMethod.q_mosaic, None, '2021-10-01'),
         ('gedi_image_list', CompositeMethod.mosaic, 'region_10000ha', None),
         ('gedi_image_list', CompositeMethod.mosaic, None, '2020-09-01'),
         ('l8_9_image_list', CompositeMethod.medoid, 'region_10000ha', None),
         ('l8_9_image_list', CompositeMethod.medoid, None, '2021-10-01'),
-    ]
+    ],
 )
 def test_composite_region_date_ordering(image_list, method, region, date, request):
     """
@@ -383,13 +390,15 @@ def test_composite_region_date_ordering(image_list, method, region, date, reques
     ee_collection = gd_collection._prepare_for_composite(method=method, date=date, region=region)
     properties = gd_collection._get_properties(ee_collection)
     assert len(properties) == len(image_list)
+
     if region:
-        # test images are ordered by CLOUDLESS/FILL_PORTION
+        # test images are ordered by CLOUDLESS/FILL_PORTION, and that portions are not the same
         schema_keys = list(gd_collection.schema.keys())
-        # CLOUDLESS_PORTION is not in MaskedCollection.schema for generic images, so use FILL_PORTION instead
         portion_key = 'CLOUDLESS_PORTION' if 'CLOUDLESS_PORTION' in schema_keys else 'FILL_PORTION'
         im_portions = [im_props[portion_key] for im_props in properties.values()]
         assert sorted(im_portions) == im_portions
+        assert len(set(im_portions)) == len(im_portions)
+
     elif date:
         # test images are ordered by time difference with `date`
         im_dates = np.array(
@@ -400,60 +409,134 @@ def test_composite_region_date_ordering(image_list, method, region, date, reques
         assert all(sorted(im_date_diff, reverse=True) == im_date_diff)
 
 
+def _get_masked_portion(
+    ee_image: ee.Image,
+    proj: ee.Projection = None,
+    region: dict = None,
+) -> ee.Number:
+    """Return the valid portion of the ``ee_image`` inside ``region``.  Assumes the ``region`` is completely covered
+    by ``ee_image``.
+    """
+    proj = proj or get_projection(ee_image)
+    ee_mask = ee_image.select('SR_B.*|B.*|rh.*').mask().reduce(ee.Reducer.allNonZero()).rename('EE_MASK')
+    mean = ee_mask.reduceRegion(reducer='mean', crs=proj, scale=proj.nominalScale(), geometry=region)
+    return ee.Number(mean.get('EE_MASK')).multiply(100)
+
+
 @pytest.mark.parametrize(
-    'image_list, method, mask', [
-        ('s2_sr_image_list', CompositeMethod.q_mosaic, True), ('s2_sr_image_list', CompositeMethod.mosaic, False),
-        ('l8_9_image_list', CompositeMethod.medoid, True), ('l8_9_image_list', CompositeMethod.median, False),
-        ('s2_sr_image_list', CompositeMethod.medoid, True), ('s2_sr_image_list', CompositeMethod.medoid, False),
+    'image_list, method, mask',
+    [
+        ('s2_sr_image_list', CompositeMethod.q_mosaic, True),
+        ('s2_sr_image_list', CompositeMethod.mosaic, False),
+        ('l8_9_image_list', CompositeMethod.medoid, True),
+        ('l8_9_image_list', CompositeMethod.median, False),
+        ('s2_sr_image_list', CompositeMethod.medoid, True),
+        ('s2_sr_image_list', CompositeMethod.medoid, False),
         ('l8_9_image_list', CompositeMethod.medoid, False),
-    ]
+    ],
 )
 def test_composite_mask(image_list, method, mask, region_100ha, request):
-    """
-    Test the MaskedImage.composite() `mask` parameter results in a masked composite image, comprised of masked
-    component images.
-    """
-
+    """In MaskedImage.composite(), test masking of component and composite images with the `mask` parameter."""
+    # form the composite collection and image
     image_list: List = request.getfixturevalue(image_list)
     gd_collection = MaskedCollection.from_list(image_list)
-    proj = get_projection(gd_collection.ee_collection.first(), min_scale=True)
     ee_collection = gd_collection._prepare_for_composite(method=method, mask=mask)
+    composite_im = gd_collection.composite(method=method, mask=mask)
     properties = gd_collection._get_properties(ee_collection)
+    proj = get_projection(ee_collection.first(), min_scale=True)
     assert len(properties) == len(image_list)
 
-    def count_masked_pixels(ee_image: ee.Image, count_list: ee.List):
-        """ Return the pixel count (area) of the EE mask / valid image area.  """
-        ee_mask = ee_image.select('SR_B.*|B.*|rh.*').mask().reduce(ee.Reducer.allNonZero()).rename('EE_MASK')
-        count = ee_mask.reduceRegion(reducer='sum', crs=proj.crs(), scale=proj.nominalScale(), geometry=region_100ha)
-        return ee.List(count_list).add(count.get('EE_MASK'))
+    def get_masked_portions(ee_image: ee.Image, portions: ee.List) -> ee.List:
+        """Append the valid portion of ``ee_image`` to  ``portions``."""
+        portion = _get_masked_portion(ee_image, proj=proj, region=region_100ha)
+        return ee.List(portions).add(portion)
 
-    # get the mask pixel counts for the component images
-    component_mask_counts = ee_collection.iterate(count_masked_pixels, ee.List([])).getInfo()
-    # get the mask pixel count for the composite image
-    composite_im = gd_collection.composite(method=method, mask=mask)
-    composite_mask_count = count_masked_pixels(composite_im.ee_image, []).get(0).getInfo()
+    # get the mask portions for the component and composite images
+    component_portions = np.array(ee_collection.iterate(get_masked_portions, ee.List([])).getInfo())
+    composite_portion = _get_masked_portion(composite_im.ee_image, proj=proj, region=region_100ha).getInfo()
+
+    # test masking of components and composite image
     if mask:
-        # test the composite mask (valid area) is smaller than combined area of the component image masks
-        assert composite_mask_count <= np.sum(component_mask_counts)
-        # test the composite mask is larger than the smallest of the component image masks
-        assert composite_mask_count >= np.min(component_mask_counts)
+        assert np.all(component_portions > 0) and np.all(component_portions < 100)
+        assert composite_portion <= np.sum(component_portions)
+        assert composite_portion >= np.min(component_portions)
     else:
-        # test that the component mask areas are ~ equal to composite mask area (i.e. the image area)
-        assert np.array(component_mask_counts) == pytest.approx(composite_mask_count, rel=.01)
+        assert component_portions == pytest.approx(100, abs=1)
+        assert composite_portion >= component_portions.max()
 
 
 @pytest.mark.parametrize(
-    'image_list, resampling, std_scale', [
-        ('s2_sr_image_list', ResamplingMethod.bilinear, 60), ('s2_sr_image_list', ResamplingMethod.bicubic, 60),
-        ('s2_sr_image_list', ResamplingMethod.average, 120), ('l8_9_image_list', ResamplingMethod.bilinear, 30),
-        ('l8_9_image_list', ResamplingMethod.bicubic, 30), ('l8_9_image_list', ResamplingMethod.average, 120),
-    ]
+    'masked_image, method, mask_kwargs',
+    [
+        ('s2_sr_hm_nocp_masked_image', 'q-mosaic', dict(mask_method='cloud-prob')),
+        ('s2_sr_hm_qa_mask_masked_image', 'mosaic', dict(mask_method='qa')),
+        ('s2_sr_hm_qa_zero_masked_image', 'medoid', dict(mask_method='qa')),
+    ],
+)
+def test_s2_composite_mask_missing_data(masked_image: str, method: str, mask_kwargs: dict, region_100ha, request):
+    """In MaskedImage.composite(), test when an S2 component image is fully masked due to missing cloud data,
+    the composite image is also fully masked.
+    """
+    # form the composite collection and image
+    masked_image: List = request.getfixturevalue(masked_image)
+    image_list = [masked_image]
+    gd_collection = MaskedCollection.from_list(image_list)
+
+    ee_collection = gd_collection._prepare_for_composite(method=method, mask=True, **mask_kwargs)
+    composite_im = gd_collection.composite(method=method, mask=True, **mask_kwargs)
+    properties = gd_collection._get_properties(ee_collection)
+    proj = get_projection(ee_collection.first(), min_scale=True)
+    assert len(properties) == len(image_list)
+
+    # get the mask portions for the component and composite images
+    component_portion = _get_masked_portion(ee_collection.first(), proj=proj, region=region_100ha).getInfo()
+    composite_portion = _get_masked_portion(composite_im.ee_image, proj=proj, region=region_100ha).getInfo()
+
+    # test component and composite images are fully masked
+    assert component_portion == composite_portion == 0
+
+
+def test_s2_composite_q_mosaic_missing_data(s2_sr_hm_nocp_masked_image: MaskedImage, region_100ha: dict):
+    """In MaskedImage.composite(), test when an S2 component image is unmasked, but has masked CLOUD_DIST band due to
+    missing cloud data, the composite image is fully masked with 'q-mosaic' method.
+    """
+    # form the composite collection and image
+    image_list = [s2_sr_hm_nocp_masked_image]
+    gd_collection = MaskedCollection.from_list(image_list)
+
+    kwargs = dict(method='q-mosaic', mask_method='cloud-prob', mask=False)
+    ee_collection = gd_collection._prepare_for_composite(**kwargs)
+    composite_im = gd_collection.composite(**kwargs)
+    properties = gd_collection._get_properties(ee_collection)
+    proj = get_projection(ee_collection.first(), min_scale=True)
+    assert len(properties) == len(image_list)
+
+    # get and test the mask portions for the component and composite images
+    component_portion = _get_masked_portion(ee_collection.first(), proj=proj, region=region_100ha).getInfo()
+    composite_portion = _get_masked_portion(composite_im.ee_image, proj=proj, region=region_100ha).getInfo()
+    assert component_portion == pytest.approx(100, abs=1)
+    assert composite_portion == 0
+
+
+@pytest.mark.parametrize(
+    'image_list, resampling, std_scale',
+    [
+        ('s2_sr_image_list', ResamplingMethod.bilinear, 60),
+        ('s2_sr_image_list', ResamplingMethod.bicubic, 60),
+        ('s2_sr_image_list', ResamplingMethod.average, 120),
+        ('l8_9_image_list', ResamplingMethod.bilinear, 30),
+        ('l8_9_image_list', ResamplingMethod.bicubic, 30),
+        ('l8_9_image_list', ResamplingMethod.average, 120),
+    ],
 )
 def test_composite_resampling(
-    image_list: str, resampling: ResamplingMethod, std_scale: float, region_10000ha: Dict,
-    request: pytest.FixtureRequest
+    image_list: str,
+    resampling: ResamplingMethod,
+    std_scale: float,
+    region_10000ha: Dict,
+    request: pytest.FixtureRequest,
 ):
-    """ Test that resampling smooths the composite image. """
+    """Test that resampling smooths the composite image."""
 
     image_list: List = request.getfixturevalue(image_list)
     gd_collection = MaskedCollection.from_list(image_list)
@@ -497,16 +580,28 @@ def test_composite_landsat_cloud_mask_params(l8_9_image_list, region_10000ha):
 
 
 @pytest.mark.parametrize(
-    'image_list, method, mask, region, date, cloud_kwargs', [
+    'image_list, method, mask, region, date, cloud_kwargs',
+    [
         ('s2_sr_image_list', CompositeMethod.q_mosaic, True, 'region_100ha', None, {}),
         ('s2_sr_image_list', CompositeMethod.mosaic, True, None, '2021-10-01', {}),
         ('s2_sr_image_list', CompositeMethod.medoid, False, None, None, {}),
         (
-            's2_sr_image_list', CompositeMethod.median, True, None, None,
+            's2_sr_image_list',
+            CompositeMethod.median,
+            True,
+            None,
+            None,
             dict(
-                mask_method='qa', mask_cirrus=False, mask_shadows=False, prob=60, dark=0.2, shadow_dist=500, buffer=500,
-                cdi_thresh=None, max_cloud_dist=2000
-            )
+                mask_method='qa',
+                mask_cirrus=False,
+                mask_shadows=False,
+                prob=60,
+                dark=0.2,
+                shadow_dist=500,
+                buffer=500,
+                cdi_thresh=None,
+                max_cloud_dist=2000,
+            ),
         ),
         ('l8_9_image_list', CompositeMethod.q_mosaic, True, 'region_100ha', None, {}),
         ('l8_9_image_list', CompositeMethod.mosaic, False, None, '2022-03-01', {}),
@@ -514,18 +609,22 @@ def test_composite_landsat_cloud_mask_params(l8_9_image_list, region_10000ha):
         ('l8_9_image_list', CompositeMethod.median, True, None, None, {}),
         ('l4_5_image_list', CompositeMethod.q_mosaic, False, 'region_100ha', None, {}),
         (
-            'l4_5_image_list', CompositeMethod.mosaic, True, None, '1988-01-01',
-            dict(mask_cirrus=False, mask_shadows=False)
+            'l4_5_image_list',
+            CompositeMethod.mosaic,
+            True,
+            None,
+            '1988-01-01',
+            dict(mask_cirrus=False, mask_shadows=False),
         ),
         ('l4_5_image_list', CompositeMethod.medoid, True, None, None, {}),
         ('l4_5_image_list', CompositeMethod.median, True, None, None, {}),
         ('gedi_image_list', CompositeMethod.mosaic, True, 'region_100ha', None, {}),
         ('gedi_image_list', CompositeMethod.mosaic, True, None, '2020-09-01', {}),
         ('gedi_image_list', CompositeMethod.medoid, True, None, None, {}),
-    ]
-)  # yapf: disable
+    ],
+)
 def test_composite(image_list, method, mask, region, date, cloud_kwargs, request):
-    """ Test MaskedCollection.composite() runs successfully with a variety of `method` and other parameters. """
+    """Test MaskedCollection.composite() runs successfully with a variety of `method` and other parameters."""
     image_list: List = request.getfixturevalue(image_list)
     region: Dict = request.getfixturevalue(region) if region else None
     gd_collection = MaskedCollection.from_list(image_list)
@@ -535,7 +634,7 @@ def test_composite(image_list, method, mask, region, date, cloud_kwargs, request
 
 
 def test_composite_errors(gedi_image_list, region_100ha):
-    """ Test MaskedCollection.composite() error conditions. """
+    """Test MaskedCollection.composite() error conditions."""
     gedi_collection = MaskedCollection.from_list(gedi_image_list)
     with pytest.raises(ValueError):
         # q-mosaic is only supported on cloud/shadow maskable images
@@ -554,14 +653,14 @@ def test_composite_errors(gedi_image_list, region_100ha):
 
 @pytest.mark.parametrize('image_list', ['s2_sr_image_list', 'l8_9_image_list'])
 def test_composite_date(image_list: str, request: pytest.FixtureRequest):
-    """ Test the composite date is the same as the first input image date. """
+    """Test the composite date is the same as the first input image date."""
 
     image_list: List = request.getfixturevalue(image_list)
     gd_collection = MaskedCollection.from_list(image_list)
     # assumes the image_list's are in date order
     first_date = datetime.utcfromtimestamp(
         gd_collection.ee_collection.first().get('system:time_start').getInfo() / 1000
-    ) # yapf: disable
+    )
     comp_im = gd_collection.composite()
     assert comp_im.date == first_date
 
@@ -585,9 +684,7 @@ def test_composite_mult_kwargs(region_100ha):
     assert cp_prob80 != pytest.approx(cp_prob40, abs=1e-1)
 
 
-@pytest.mark.parametrize(
-    'name', ['FAO/WAPOR/2/L1_RET_E', 'MODIS/006/MCD43A4']
-)
+@pytest.mark.parametrize('name', ['FAO/WAPOR/2/L1_RET_E', 'MODIS/006/MCD43A4'])
 def test_unbounded_search_no_region(name):
     """
     Test searching an unbounded collection without a region raises an exception.
@@ -612,4 +709,3 @@ def test_unknown_collection_error(region_25ha):
         gd_collection = gd_collection.search(start_date, end_date, region_25ha)
         _ = gd_collection.properties
     assert 'not found' in str(ex)
-
diff --git a/tests/test_download.py b/tests/test_download.py
index 6516212..d01dd42 100644
--- a/tests/test_download.py
+++ b/tests/test_download.py
@@ -127,7 +127,6 @@ def _test_export_profile(exp_profile, tgt_profile, transform_shape=False):
         assert exp_profile[key] == tgt_profile[key]
 
     assert utils.rio_crs(exp_profile['crs']) == utils.rio_crs(tgt_profile['crs'])
-
     assert exp_profile['transform'][0] == tgt_profile['transform'][0]
 
     if transform_shape:
@@ -466,18 +465,15 @@ def test_prepare_for_download(base_image: str, request: pytest.FixtureRequest):
         ('int16', -(2**15)),
         ('uint32', 0),
         ('int32', -(2**31)),
-        ('float32', float('nan')),
-        ('float64', float('nan')),
+        ('float32', float('-inf')),
+        ('float64', float('-inf')),
     ],
 )
 def test_prepare_nodata(user_fix_base_image: BaseImage, region_25ha: Dict, dtype: str, exp_nodata: float):
     """Test BaseImage._prepare_for_download() sets rasterio profile nodata correctly for different dtypes."""
     exp_image, exp_profile = user_fix_base_image._prepare_for_download(region=region_25ha, scale=30, dtype=dtype)
     assert exp_image.dtype == dtype
-    if np.isnan(exp_profile['nodata']):
-        assert np.isnan(exp_nodata)
-    else:
-        assert exp_profile['nodata'] == exp_nodata
+    assert exp_profile['nodata'] == exp_nodata
 
 
 @pytest.mark.parametrize(
diff --git a/tests/test_mask.py b/tests/test_mask.py
index 84d30fd..6f530ef 100644
--- a/tests/test_mask.py
+++ b/tests/test_mask.py
@@ -13,17 +13,20 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
-from typing import Dict
+
+from typing import Dict, Iterable
 
 import ee
 import numpy as np
 import pytest
 import rasterio as rio
-from geedim.mask import MaskedImage, get_projection, class_from_id
+
+from geedim.enums import CloudMaskMethod
+from geedim.mask import class_from_id, CloudMaskedImage, MaskedImage, Sentinel2SrClImage
 
 
 def test_class_from_id(landsat_image_ids, s2_sr_image_id, s2_toa_hm_image_id, generic_image_ids):
-    """  Test class_from_id(). """
+    """Test class_from_id()."""
     from geedim.mask import LandsatImage, Sentinel2SrClImage, Sentinel2ToaClImage
 
     assert all([class_from_id(im_id) == LandsatImage for im_id in landsat_image_ids])
@@ -33,33 +36,47 @@ def test_class_from_id(landsat_image_ids, s2_sr_image_id, s2_toa_hm_image_id, ge
 
 
 def test_from_id():
-    """ Test MaskedImage.from_id() sets _id. """
+    """Test MaskedImage.from_id() sets _id."""
     ee_id = 'MODIS/006/MCD43A4/2022_01_01'
     gd_image = MaskedImage.from_id(ee_id)
     assert gd_image._id == ee_id
 
 
 @pytest.mark.parametrize(
-    'masked_image', [
-        'user_masked_image', 'modis_nbar_masked_image', 'gch_masked_image', 's1_sar_masked_image',
-        'gedi_agb_masked_image', 'gedi_cth_masked_image', 'landsat_ndvi_masked_image'
-    ]
+    'masked_image',
+    [
+        'user_masked_image',
+        'modis_nbar_masked_image',
+        'gch_masked_image',
+        's1_sar_masked_image',
+        'gedi_agb_masked_image',
+        'gedi_cth_masked_image',
+        'landsat_ndvi_masked_image',
+    ],
 )
 def test_gen_aux_bands_exist(masked_image: str, request: pytest.FixtureRequest):
-    """ Test the presence of auxiliary band (i.e. FILL_MASK) in generic masked images. """
+    """Test the presence of auxiliary band (i.e. FILL_MASK) in generic masked images."""
     masked_image: MaskedImage = request.getfixturevalue(masked_image)
     band_names = masked_image.ee_image.bandNames().getInfo()
     assert 'FILL_MASK' in band_names
 
 
 @pytest.mark.parametrize(
-    'masked_image', [
-        's2_sr_masked_image', 's2_toa_masked_image', 's2_sr_hm_masked_image', 's2_toa_hm_masked_image',
-        'l9_masked_image', 'l8_masked_image', 'l7_masked_image', 'l5_masked_image', 'l4_masked_image'
-    ]
+    'masked_image',
+    [
+        's2_sr_masked_image',
+        's2_toa_masked_image',
+        's2_sr_hm_masked_image',
+        's2_toa_hm_masked_image',
+        'l9_masked_image',
+        'l8_masked_image',
+        'l7_masked_image',
+        'l5_masked_image',
+        'l4_masked_image',
+    ],
 )
 def test_cloud_mask_aux_bands_exist(masked_image: str, request: pytest.FixtureRequest):
-    """ Test the presence of auxiliary bands in cloud masked images. """
+    """Test the presence of auxiliary bands in cloud masked images."""
     masked_image: MaskedImage = request.getfixturevalue(masked_image)
     band_names = masked_image.ee_image.bandNames().getInfo()
     exp_band_names = ['CLOUD_MASK', 'SHADOW_MASK', 'FILL_MASK', 'CLOUDLESS_MASK', 'CLOUD_DIST']
@@ -68,27 +85,63 @@ def test_cloud_mask_aux_bands_exist(masked_image: str, request: pytest.FixtureRe
 
 
 @pytest.mark.parametrize(
-    'masked_image', [
-        's2_sr_masked_image', 's2_toa_masked_image', 's2_sr_hm_masked_image', 's2_toa_hm_masked_image',
-        'l9_masked_image', 'l8_masked_image', 'l7_masked_image', 'l5_masked_image', 'l4_masked_image',
-        'user_masked_image', 'modis_nbar_masked_image', 'gch_masked_image', 's1_sar_masked_image',
-        'gedi_agb_masked_image', 'gedi_cth_masked_image', 'landsat_ndvi_masked_image'
-    ]
+    'masked_image',
+    [
+        's2_sr_masked_image',
+        's2_toa_masked_image',
+        's2_sr_hm_masked_image',
+        's2_toa_hm_masked_image',
+        'l9_masked_image',
+        'l8_masked_image',
+        'l7_masked_image',
+        'l5_masked_image',
+        'l4_masked_image',
+        'user_masked_image',
+        'modis_nbar_masked_image',
+        'gch_masked_image',
+        's1_sar_masked_image',
+        'gedi_agb_masked_image',
+        'gedi_cth_masked_image',
+        'landsat_ndvi_masked_image',
+    ],
 )
 def test_set_region_stats(masked_image: str, region_100ha, request: pytest.FixtureRequest):
     """
     Test MaskedImage._set_region_stats() generates the expected properties and that these are in the valid range.
     """
     masked_image: MaskedImage = request.getfixturevalue(masked_image)
-    masked_image._set_region_stats(region_100ha)
+    masked_image._set_region_stats(region_100ha, scale=masked_image._ee_proj.nominalScale())
     for stat_name in ['FILL_PORTION', 'CLOUDLESS_PORTION']:
         assert stat_name in masked_image.properties
         assert masked_image.properties[stat_name] >= 0 and masked_image.properties[stat_name] <= 100
 
 
+@pytest.mark.parametrize(
+    'masked_image, exp_scale',
+    [
+        ('s2_sr_hm_masked_image', 60),
+        ('l9_masked_image', 30),
+        ('l4_masked_image', 30),
+        ('s1_sar_masked_image', 10),
+        ('gedi_agb_masked_image', 1000),
+        # include fixtures with bands that have no fixed projection
+        ('s2_sr_hm_qa_mask_masked_image', 60),
+        ('s2_sr_hm_qa_zero_masked_image', 60),
+        ('s2_sr_hm_nocp_masked_image', 60),
+    ],
+)
+def test_ee_proj(masked_image: str, exp_scale: float, request: pytest.FixtureRequest):
+    """Test MaskedImage._ee_proj has the correct scale and CRS."""
+    masked_image: MaskedImage = request.getfixturevalue(masked_image)
+    proj = masked_image._ee_proj.getInfo()
+
+    assert np.abs(proj['transform'][0]) == pytest.approx(exp_scale, rel=1e-3)
+    assert proj.get('crs', 'wkt') != 'EPSG:4326'
+
+
 @pytest.mark.parametrize('image_id', ['l9_image_id', 'l8_image_id', 'l7_image_id', 'l5_image_id', 'l4_image_id'])
 def test_landsat_cloudless_portion(image_id: str, request: pytest.FixtureRequest):
-    """ Test `geedim` CLOUDLESS_PORTION for the whole image against related Landsat CLOUD_COVER property. """
+    """Test `geedim` CLOUDLESS_PORTION for the whole image against related Landsat CLOUD_COVER property."""
     image_id: MaskedImage = request.getfixturevalue(image_id)
     masked_image = MaskedImage.from_id(image_id, mask_shadows=False, mask_cirrus=False)
     masked_image._set_region_stats()
@@ -101,7 +154,7 @@ def test_landsat_cloudless_portion(image_id: str, request: pytest.FixtureRequest
 
 @pytest.mark.parametrize('image_id', ['s2_toa_image_id', 's2_sr_image_id', 's2_toa_hm_image_id', 's2_sr_hm_image_id'])
 def test_s2_cloudless_portion(image_id: str, request: pytest.FixtureRequest):
-    """ Test `geedim` CLOUDLESS_PORTION for the whole image against CLOUDY_PIXEL_PERCENTAGE Sentinel-2 property. """
+    """Test `geedim` CLOUDLESS_PORTION for the whole image against CLOUDY_PIXEL_PERCENTAGE Sentinel-2 property."""
     image_id: str = request.getfixturevalue(image_id)
     masked_image = MaskedImage.from_id(image_id, mask_method='qa', mask_shadows=False, mask_cirrus=False)
     masked_image._set_region_stats()
@@ -113,7 +166,7 @@ def test_s2_cloudless_portion(image_id: str, request: pytest.FixtureRequest):
 
 @pytest.mark.parametrize('image_id', ['l9_image_id'])
 def test_landsat_cloudmask_params(image_id: str, request: pytest.FixtureRequest):
-    """ Test Landsat cloud/shadow masking `mask_shadows` and `mask_cirrus` parameters. """
+    """Test Landsat cloud/shadow masking `mask_shadows` and `mask_cirrus` parameters."""
     image_id: str = request.getfixturevalue(image_id)
     masked_image = MaskedImage.from_id(image_id, mask_shadows=False, mask_cirrus=False)
     masked_image._set_region_stats()
@@ -135,7 +188,7 @@ def test_landsat_cloudmask_params(image_id: str, request: pytest.FixtureRequest)
 
 @pytest.mark.parametrize('image_id', ['s2_sr_image_id', 's2_toa_image_id'])
 def test_s2_cloudmask_mask_shadows(image_id: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test S2 cloud/shadow masking `mask_shadows` parameter. """
+    """Test S2 cloud/shadow masking `mask_shadows` parameter."""
     image_id: str = request.getfixturevalue(image_id)
     masked_image = MaskedImage.from_id(image_id, mask_shadows=False)
     masked_image._set_region_stats(region_10000ha)
@@ -150,7 +203,7 @@ def test_s2_cloudmask_mask_shadows(image_id: str, region_10000ha: Dict, request:
 
 @pytest.mark.parametrize('image_id', ['s2_sr_image_id', 's2_toa_image_id'])
 def test_s2_cloudmask_prob(image_id: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test S2 cloud/shadow masking `prob` parameter. """
+    """Test S2 cloud/shadow masking `prob` parameter."""
     image_id: str = request.getfixturevalue(image_id)
     masked_image = MaskedImage.from_id(image_id, mask_shadows=True, prob=80)
     masked_image._set_region_stats(region_10000ha)
@@ -159,27 +212,28 @@ def test_s2_cloudmask_prob(image_id: str, region_10000ha: Dict, request: pytest.
     masked_image._set_region_stats(region_10000ha)
     prob40_portion = 100 * masked_image.properties['CLOUDLESS_PORTION'] / masked_image.properties['FILL_PORTION']
     # test there is more cloud (less CLOUDLESS_PORTION) with prob=40 as compared to prob=80
-    assert prob80_portion > prob40_portion
+    assert prob80_portion > prob40_portion > 0
 
 
 @pytest.mark.parametrize('image_id', ['s2_sr_image_id', 's2_toa_image_id'])
 def test_s2_cloudmask_method(image_id: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test S2 cloud/shadow masking `mask_method` parameter. """
+    """Test S2 cloud/shadow masking `mask_method` parameter."""
     image_id: str = request.getfixturevalue(image_id)
-    masked_image = MaskedImage.from_id(image_id, mask_method='cloud-prob')
-    masked_image._set_region_stats(region_10000ha)
-    cloud_prob_portion = 100 * masked_image.properties['CLOUDLESS_PORTION'] / masked_image.properties['FILL_PORTION']
-    masked_image = MaskedImage.from_id(image_id, mask_method='qa')
-    masked_image._set_region_stats(region_10000ha)
-    qa_portion = 100 * masked_image.properties['CLOUDLESS_PORTION'] / masked_image.properties['FILL_PORTION']
+    cl_portions = []
+    for mask_method in CloudMaskMethod:
+        masked_image = MaskedImage.from_id(image_id, mask_method=mask_method)
+        masked_image._set_region_stats(region_10000ha)
+        cl_portions.append(100 * masked_image.properties['CLOUDLESS_PORTION'])
+
     # test `mask_method` changes CLOUDLESS_PORTION but not by too much
-    assert cloud_prob_portion != qa_portion
-    assert cloud_prob_portion == pytest.approx(qa_portion, abs=10)
+    assert len(set(cl_portions)) == len(cl_portions)
+    assert all([cl_portions[0] != pytest.approx(cp, abs=10) for cp in cl_portions[1:]])
+    assert all([cp != pytest.approx(0, abs=1) for cp in cl_portions])
 
 
 @pytest.mark.parametrize('image_id', ['s2_sr_image_id', 's2_toa_image_id'])
 def test_s2_cloudmask_mask_cirrus(image_id: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test S2 cloud/shadow masking `mask_cirrus` parameter. """
+    """Test S2 cloud/shadow masking `mask_cirrus` parameter."""
     image_id: str = request.getfixturevalue(image_id)
     masked_image = MaskedImage.from_id(image_id, mask_method='qa', mask_cirrus=False)
     # cloud and shadow free portion
@@ -194,7 +248,7 @@ def test_s2_cloudmask_mask_cirrus(image_id: str, region_10000ha: Dict, request:
 
 @pytest.mark.parametrize('image_id', ['s2_sr_image_id', 's2_toa_image_id'])
 def test_s2_cloudmask_dark(image_id: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test S2 cloud/shadow masking `dark` parameter. """
+    """Test S2 cloud/shadow masking `dark` parameter."""
     image_id: str = request.getfixturevalue(image_id)
     masked_image = MaskedImage.from_id(image_id, dark=0.5)
     masked_image._set_region_stats(region_10000ha)
@@ -209,7 +263,7 @@ def test_s2_cloudmask_dark(image_id: str, region_10000ha: Dict, request: pytest.
 
 @pytest.mark.parametrize('image_id', ['s2_sr_image_id', 's2_toa_image_id'])
 def test_s2_cloudmask_shadow_dist(image_id: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test S2 cloud/shadow masking `shadow_dist` parameter. """
+    """Test S2 cloud/shadow masking `shadow_dist` parameter."""
     image_id: str = request.getfixturevalue(image_id)
     masked_image = MaskedImage.from_id(image_id, shadow_dist=200)
     masked_image._set_region_stats(region_10000ha)
@@ -224,12 +278,12 @@ def test_s2_cloudmask_shadow_dist(image_id: str, region_10000ha: Dict, request:
 
 @pytest.mark.parametrize('image_id', ['s2_sr_image_id', 's2_toa_image_id'])
 def test_s2_cloudmask_cdi_thresh(image_id: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test S2 cloud/shadow masking `cdi_thresh` parameter. """
+    """Test S2 cloud/shadow masking `cdi_thresh` parameter."""
     image_id: str = request.getfixturevalue(image_id)
-    masked_image = MaskedImage.from_id(image_id, cdi_thresh=.5)
+    masked_image = MaskedImage.from_id(image_id, cdi_thresh=0.5)
     masked_image._set_region_stats(region_10000ha)
     cdi_pt5_portion = 100 * masked_image.properties['CLOUDLESS_PORTION'] / masked_image.properties['FILL_PORTION']
-    masked_image = MaskedImage.from_id(image_id, cdi_thresh=-.5)
+    masked_image = MaskedImage.from_id(image_id, cdi_thresh=-0.5)
     masked_image._set_region_stats(region_10000ha)
     cdi_negpt5_portion = 100 * masked_image.properties['CLOUDLESS_PORTION'] / masked_image.properties['FILL_PORTION']
     # test that increasing `cdi_thresh` results in an increase in detected cloud and corresponding decrease in
@@ -239,10 +293,10 @@ def test_s2_cloudmask_cdi_thresh(image_id: str, region_10000ha: Dict, request: p
 
 @pytest.mark.parametrize('image_id, max_cloud_dist', [('s2_sr_image_id', 100), ('s2_sr_hm_image_id', 500)])
 def test_s2_clouddist_max(image_id: str, max_cloud_dist: int, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test S2 cloud distance `max_cloud_dist` parameter. """
+    """Test S2 cloud distance `max_cloud_dist` parameter."""
 
     def get_max_cloud_dist(cloud_dist: ee.Image):
-        """ Get the maximum of `cloud_dist` over region_10000ha. """
+        """Get the maximum of `cloud_dist` over region_10000ha."""
         mcd = cloud_dist.reduceRegion(reducer='max', geometry=region_10000ha, bestEffort=True, maxPixels=1e4)
         return mcd.get('CLOUD_DIST').getInfo() * 10
 
@@ -253,82 +307,170 @@ def get_max_cloud_dist(cloud_dist: ee.Image):
     assert meas_max_cloud_dist == pytest.approx(max_cloud_dist, rel=0.1)
 
 
-def get_mask_stats(masked_image: MaskedImage, region: Dict):
-    """  Get cloud/shadow etc aux band statistics for a given MaskedImage and region. """
+@pytest.mark.parametrize(
+    'masked_image', ['s2_sr_hm_qa_mask_masked_image', 's2_sr_hm_qa_zero_masked_image', 's2_sr_hm_nocp_masked_image']
+)
+def test_s2_region_stats_missing_data(masked_image: str, region_10000ha: dict, request: pytest.FixtureRequest):
+    """Test S2 region stats for unmasked images missing required cloud data."""
+    masked_image: MaskedImage = request.getfixturevalue(masked_image)
+    masked_image._set_region_stats(region_10000ha, scale=60)
+
+    assert masked_image.properties is not None
+    assert masked_image.properties['CLOUDLESS_PORTION'] == pytest.approx(0, abs=1)
+    assert masked_image.properties['FILL_PORTION'] == pytest.approx(100, abs=1)
+
+
+def _test_aux_stats(masked_image: MaskedImage, region: Dict):
+    """Sanity tests on cloud/shadow etc. aux bands for a given MaskedImage and region."""
+    # create a pan image for testing brightnesses
     ee_image = masked_image.ee_image
     pan = ee_image.select([0, 1, 2]).reduce(ee.Reducer.mean())
-    cdist = ee_image.select('CLOUD_DIST')
+
+    # mask the pan image with cloud, shadow & cloudless masks
     cloud_mask = ee_image.select('CLOUD_MASK')
     shadow_mask = ee_image.select('SHADOW_MASK')
     cloudless_mask = ee_image.select('CLOUDLESS_MASK')
     pan_cloud = pan.updateMask(cloud_mask).rename('PAN_CLOUD')
-    pan_shadow = pan.mask(shadow_mask).rename('PAN_SHADOW')
+    pan_shadow = pan.updateMask(shadow_mask).rename('PAN_SHADOW')
     pan_cloudless = pan.updateMask(cloudless_mask).rename('PAN_CLOUDLESS')
+
+    # mask the cloud distance image with cloud & cloudless masks
+    cdist = ee_image.select('CLOUD_DIST')
     cdist_cloud = cdist.updateMask(cloud_mask).rename('CDIST_CLOUD')
     cdist_cloudless = cdist.updateMask(cloudless_mask).rename('CDIST_CLOUDLESS')
+
+    # find mean stats of all masked images, and min of cloud distance where it is >0
     stats_image = ee.Image([pan_cloud, pan_shadow, pan_cloudless, cdist_cloud, cdist_cloudless])
-    proj = get_projection(ee_image, min_scale=False)
-    means = stats_image.reduceRegion(
-        reducer='mean', geometry=region, crs=proj.crs(), scale=proj.nominalScale(), bestEffort=True, maxPixels=1e8
+    proj = masked_image._ee_proj
+    stats = stats_image.reduceRegion(
+        reducer='mean', geometry=region, crs=proj, scale=proj.nominalScale(), bestEffort=True, maxPixels=1e8
     )
     cdist_min = cdist.updateMask(cdist).reduceRegion(
-        reducer='min', geometry=region, crs=proj.crs(), scale=proj.nominalScale(), bestEffort=True, maxPixels=1e8
+        reducer='min', geometry=region, crs=proj, scale=proj.nominalScale(), bestEffort=True, maxPixels=1e8
     )
-    means = means.set('CDIST_MIN', cdist_min.get('CLOUD_DIST'))
-    return means.getInfo()
+    stats = stats.set('CDIST_MIN', cdist_min.get('CLOUD_DIST'))
+    stats = stats.getInfo()
+
+    # test cloud is brighter than cloudless
+    assert stats['PAN_CLOUD'] > stats['PAN_CLOUDLESS']
+    # test cloudless is brighter than shadow
+    assert stats['PAN_CLOUDLESS'] > stats['PAN_SHADOW']
+    # test cloudless areas have greater distance to cloud than cloudy areas
+    assert stats['CDIST_CLOUDLESS'] > stats['CDIST_CLOUD']
+    # test min distance to cloud is pixel size
+    assert stats['CDIST_MIN'] * 10 == masked_image._ee_proj.nominalScale().getInfo()
 
 
 @pytest.mark.parametrize(
     'masked_image', ['l9_masked_image', 'l8_masked_image', 'l7_masked_image', 'l5_masked_image', 'l4_masked_image']
 )
 def test_landsat_aux_bands(masked_image: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test Landsat auxiliary band values. """
+    """Test Landsat auxiliary band values for sanity."""
     masked_image: MaskedImage = request.getfixturevalue(masked_image)
-    stats = get_mask_stats(masked_image, region_10000ha)
-    assert stats['PAN_CLOUD'] > stats['PAN_CLOUDLESS']
-    assert stats['PAN_CLOUDLESS'] > stats['PAN_SHADOW']
-    assert stats['CDIST_CLOUDLESS'] > stats['CDIST_CLOUD']
+    _test_aux_stats(masked_image, region_10000ha)
 
 
-@pytest.mark.parametrize('masked_image',
-    ['s2_sr_masked_image', 's2_toa_masked_image', 's2_sr_hm_masked_image', 's2_toa_hm_masked_image']
-)  # yapf: disable
-def test_s2_aux_bands(masked_image: str, region_10000ha: Dict, request: pytest.FixtureRequest):
-    """ Test Sentinel-2 auxiliary band values. """
+@pytest.mark.parametrize(
+    'image_id, mask_methods',
+    [
+        ('s2_sr_image_id', CloudMaskMethod),
+        ('s2_toa_image_id', CloudMaskMethod),
+        ('s2_sr_hm_image_id', CloudMaskMethod),
+        ('s2_toa_hm_image_id', CloudMaskMethod),
+        # images missing QA60 so do cloud-prob method only
+        ('s2_sr_hm_qa_mask_image_id', ['cloud-prob']),
+        ('s2_sr_hm_qa_zero_image_id', ['cloud-prob']),
+    ],
+)
+def test_s2_aux_bands(image_id: str, mask_methods: Iterable, region_10000ha: Dict, request: pytest.FixtureRequest):
+    """Test Sentinel-2 auxiliary band values for sanity with all masking methods."""
+    image_id: str = request.getfixturevalue(image_id)
+    for mask_method in mask_methods:
+        masked_image = MaskedImage.from_id(image_id, mask_method=mask_method)
+        _test_aux_stats(masked_image, region_10000ha)
+
+
+def test_s2_aux_bands_unlink(s2_sr_hm_image_id: str, region_10000ha: Dict):
+    """Test Sentinel-2 auxiliary band values for sanity on an image without linked cloud data."""
+    # TODO: include the cloud score+ method in this test when that method is added
+    # create an image with unknown id to prevent linking to cloud data
+    ee_image = ee.Image(s2_sr_hm_image_id)
+    ee_image = ee_image.set('system:index', 'COPERNICUS/S2_HARMONIZED/unknown')
+
+    for mask_method in ['qa']:
+        masked_image = Sentinel2SrClImage(ee_image, mask_method=mask_method)
+        _test_aux_stats(masked_image, region_10000ha)
+
+
+@pytest.mark.parametrize(
+    'masked_image', ['s2_sr_hm_nocp_masked_image', 's2_sr_hm_qa_mask_masked_image', 's2_sr_hm_qa_zero_masked_image']
+)
+def test_s2_aux_bands_missing_data(masked_image: str, region_10000ha: Dict, request: pytest.FixtureRequest):
+    """Test Sentinel-2 auxiliary band masking / transparency for unmasked images missing required cloud data."""
     masked_image: MaskedImage = request.getfixturevalue(masked_image)
-    stats = get_mask_stats(masked_image, region_10000ha)
-    assert stats['PAN_CLOUD'] > stats['PAN_CLOUDLESS']
-    assert stats['PAN_CLOUDLESS'] > stats['PAN_SHADOW']
-    assert stats['CDIST_CLOUDLESS'] > stats['CDIST_CLOUD']
-    proj_scale = get_projection(masked_image.ee_image, min_scale=False).nominalScale().getInfo()
-    assert stats['CDIST_MIN'] * 10 == proj_scale
 
+    # get region sums of the auxiliary masks
+    proj = masked_image._ee_proj
+    aux_bands = masked_image.ee_image.select('.*MASK|CLOUD_DIST')
+    aux_mask = aux_bands.mask()
+    stats = aux_mask.reduceRegion(
+        reducer='sum', geometry=region_10000ha, crs=proj, scale=proj.nominalScale(), bestEffort=True, maxPixels=1e8
+    )
+    stats = stats.getInfo()
 
-@pytest.mark.parametrize('masked_image', ['s2_sr_masked_image', 'l9_masked_image'])
+    # test auxiliary masks are transparent
+    assert stats['FILL_MASK'] > 0
+    for band_name in ['CLOUDLESS_MASK', 'CLOUD_MASK', 'SHADOW_MASK', 'CLOUD_DIST']:
+        assert stats[band_name] == 0
+
+
+@pytest.mark.parametrize('masked_image', ['gedi_cth_masked_image', 's2_sr_masked_image', 'l9_masked_image'])
 def test_mask_clouds(masked_image: str, region_100ha: Dict, tmp_path, request: pytest.FixtureRequest):
-    """ Test MaskedImage.mask_clouds() by downloading and examining dataset masks. """
+    """Test MaskedImage.mask_clouds() masks the fill or cloudless portion by downloading and examining dataset masks."""
+    masked_image: MaskedImage = request.getfixturevalue(masked_image)
+    filename = tmp_path.joinpath(f'test_image.tif')
+    masked_image.mask_clouds()
+    proj_scale = masked_image._ee_proj.nominalScale()
+    masked_image.download(filename, region=region_100ha, dtype='float32', scale=proj_scale)
+    assert filename.exists()
+
+    with rio.open(filename, 'r') as ds:
+        mask_name = 'CLOUDLESS_MASK' if isinstance(masked_image, CloudMaskedImage) else 'FILL_MASK'
+        mask = ds.read(ds.descriptions.index(mask_name) + 1, masked=True)
+
+        # test areas outside CLOUDLESS_MASK / FILL_MASK are masked
+        assert np.all(mask == 1)
+
+        # test CLOUDLESS_MASK / FILL_MASK matches the nodata mask for each band
+        mask = mask.filled(0).astype('bool')
+        ds_masks = ds.read_masks().astype('bool')
+        assert np.all(mask == ds_masks)
+
+
+@pytest.mark.parametrize(
+    'masked_image', ['s2_sr_hm_nocp_masked_image', 's2_sr_hm_qa_mask_masked_image', 's2_sr_hm_qa_zero_masked_image']
+)
+def test_s2_mask_clouds_missing_data(masked_image: str, region_100ha: Dict, tmp_path, request: pytest.FixtureRequest):
+    """Test Sentinel2SrClImage.mask_clouds() masks the entire image when it is missing required cloud data. Downloads
+    and examines dataset masks.
+    """
     masked_image: MaskedImage = request.getfixturevalue(masked_image)
     filename = tmp_path.joinpath(f'test_image.tif')
     masked_image.mask_clouds()
-    proj_scale = get_projection(masked_image.ee_image, min_scale=False).nominalScale()
-    masked_image.download(filename, region=region_100ha, dtype='int32', scale=proj_scale)
+    proj_scale = masked_image._ee_proj.nominalScale()
+    masked_image.download(filename, region=region_100ha, dtype='float32', scale=proj_scale)
     assert filename.exists()
+
+    # test all data is masked / nodata
     with rio.open(filename, 'r') as ds:
-        ds: rio.DatasetReader = ds
-        cloudless_mask = ds.read(ds.descriptions.index('CLOUDLESS_MASK') + 1, masked=True)
-        assert np.all(cloudless_mask == 1)  # all cloud/shadow areas should be masked (0)
-        cloudless_mask = cloudless_mask.filled(0).astype('bool')  # fill nodata with 0 and cast to bool
-        # test that cloudless_mask is the same as the nodata/dataset mask for each bands
         ds_masks = ds.read_masks().astype('bool')
-        assert np.all(cloudless_mask == ds_masks)
+        assert not np.any(ds_masks)
 
 
 def test_skysat_region_stats():
-    """ Test _set_region_stats() works on SKYSAT image with no region. """
+    """Test _set_region_stats() works on SKYSAT image with no region."""
     ee_image = ee.Image('SKYSAT/GEN-A/PUBLIC/ORTHO/RGB/s02_20141004T074858Z')
     masked_image = MaskedImage(ee_image)
     masked_image._set_region_stats()
     assert 'FILL_PORTION' in masked_image.properties
     assert masked_image.properties['FILL_PORTION'] > 0.8
-##
-
diff --git a/tests/test_tile.py b/tests/test_tile.py
index 7739c96..3fb49d7 100644
--- a/tests/test_tile.py
+++ b/tests/test_tile.py
@@ -74,7 +74,7 @@ def zipped_gtiff_bytes(mock_base_image: BaseImageLike) -> bytes:
             **rio.default_gtiff_profile,
             width=mock_base_image.shape[1],
             height=mock_base_image.shape[0],
-            count=mock_base_image.count
+            count=mock_base_image.count,
         ) as ds:
             ds.write(array)
         with zipfile.ZipFile(zip_buffer, 'a', zipfile.ZIP_DEFLATED, False) as zf: