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

Issue:
Shift of data in X and Y direction by 1 pixel while downloading NASADEM data in EPSG:4326
The download function at the following link does not consider crs transformation parameter, which is causing shift in downloaded data.
https://github.com/dugalh/geedim/blob/main/geedim/download.py#L686


Discussion / hints:

Fix patch:
diff --git a/geedim/cli.py b/geedim/cli.py
index b6a47f4..73404b5 100644
--- a/geedim/cli.py
+++ b/geedim/cli.py
@@ -23,6 +23,7 @@
 from typing import List
 
 import click
+import rasterio as rio
 import rasterio.crs as rio_crs
 from click.core import ParameterSource
 from geedim import schema, Initialize, version
@@ -56,15 +57,16 @@ def get_help(self, ctx):
         """ Strip some RST markup from the help text for CLI display.  Assumes no grid tables. """
         if not hasattr(self, 'click_wrap_text'):
             self.click_wrap_text = click.formatting.wrap_text
-        # TODO: copy homonim sub_strings
         sub_strings = {
-            '\b\n': '\n\b',  # convert from RST friendly to click literal (unwrapped) block marker
-            ':option:': '',  # strip ':option:'
-            '\| ': '',  # strip RST literal (unwrapped) marker in e.g. tables and bullet lists
-            '\n\.\. _.*:\n': '',  # strip RST ref directive '\n.. <name>:\n'
-            '`(.*)<(.*)>`_': '\g<1>',  # convert from RST cross-ref '`<name> <<link>>`_' to 'name'
-            '::': ':'  # convert from RST '::' to ':'
-        }
+            '\b\n': '\n\b',                 # convert from RST friendly to click literal (unwrapped) block marker
+            r'\| ': '',                     # strip RST literal (unwrapped) marker in e.g. tables and bullet lists
+            '\n\.\. _.*:\n': '',            # strip RST ref directive '\n.. _<name>:\n'
+            '`(.*?) <(.*?)>`_': r'\g<1>',   # convert from RST cross-ref '`<name> <<link>>`_' to 'name'
+            '::': ':',                      # convert from RST '::' to ':'
+            '``(.*?)``': r'\g<1>',          # convert from RST '``literal``' to 'literal'
+            ':option:`(.*?)( <.*?>)?`': r'\g<1>',  # convert ':option:`--name <group-command --name>`' to '--name'
+            ':option:`(.*?)`': r'\g<1>',    # convert ':option:`--name`' to '--name'
+        }  # yapf: disable
 
         def reformat_text(text, width, **kwargs):
             for sub_key, sub_value in sub_strings.items():
@@ -91,6 +93,13 @@ def invoke(self, ctx):
             # append any image id's to the image_list
             ctx.obj.image_list += list(ctx.params['image_id'])
 
+        if ('like' in ctx.params) and (ctx.params['like'] is not None):
+            # populate crs, crs_transform & shape parameters from a template raster
+            with rio.open(ctx.params['like'], 'r') as im:
+                ctx.params['crs'] = f'EPSG:{im.crs.to_epsg()}'  # TODO: WKT?
+                ctx.params['crs_transform'] = im.transform
+                ctx.params['shape'] = im.shape
+
         return click.Command.invoke(self, ctx)
 
 
@@ -189,15 +198,14 @@ def _prepare_image_list(obj: SimpleNamespace, mask=False) -> List[MaskedImage, ]
             raise ValueError(f'Unsupported image object type: {type(im_obj)}')
         image_list.append(im_obj)
 
-    if obj.region is None and any([not im.has_fixed_projection for im in image_list]):
-        raise click.BadOptionUsage('region', 'One of --region or --bbox is required for a composite image.')
     return image_list
 
 
 # Define click options that are common to more than one command
 bbox_option = click.option(
     '-b', '--bbox', type=click.FLOAT, nargs=4, default=None, callback=_bbox_cb,
-    help='Region defined by WGS84 bounding box co-ordinates (left, bottom, right, top).'
+    metavar='LEFT BOTTOM RIGHT TOP',
+    help='Region defined by WGS84 bounding box co-ordinates.'
 )
 region_option = click.option(
     '-r', '--region', type=click.Path(exists=True, dir_okay=False, allow_dash=True), default=None, callback=_region_cb,
@@ -230,6 +238,20 @@ def _prepare_image_list(obj: SimpleNamespace, mask=False) -> List[MaskedImage, ]
     '-so/-nso', '--scale-offset/--no-scale-offset', default=False, show_default=True,
     help='Whether to apply any EE band scales and offsets to the image.'
 )
+crs_transform_option = click.option(
+    '-ct', '--crs-transform', type=click.FLOAT, nargs=6, default=None,
+    metavar='XSCALE XSHEAR XTRANSLATION YSHEAR YSCALE YTRANSLATION',
+    help='Six element affine transform in the download CRS.  Use with ``--shape`` to specify image '
+         'bounds and resolution.'
+)  # yapf: disable
+shape_option = click.option(
+    '-sh', '--shape', type=click.INT, nargs=2, default=None, metavar='HEIGHT WIDTH',
+    help='Image height & width dimensions (pixels).'
+)
+like_option = click.option(
+    '-l', '--like', type=click.Path(exists=True, dir_okay=False), default=None,
+    help='Template raster from which to derive ``--crs``, ``--crs-transform`` & ``--shape``.'
+)
 
 
 # geedim CLI and chained command group
@@ -246,14 +268,17 @@ def cli(ctx, verbose, quiet):
 
 
 # TODO: add clear docs on what is piped out of or into each command.
-# TODO: use cloup and linear help layout.  move lists of option values from command docting to corresponding option help
+# TODO: use cloup and linear help layout & move lists of option values from command docsting to corresponding option
+#  help
+# TODO: add RST option markup like in homonim e.g. :option:`--mask-method`, and
+#  :option:`--param-image <homonim-fuse --param-image>`
 
 # config command
 @click.command(cls=ChainedCommand, context_settings=dict(auto_envvar_prefix='GEEDIM'))
 @click.option(
     '-mc/-nmc', '--mask-cirrus/--no-mask-cirrus', default=True, show_default=True,
     help='Whether to mask cirrus clouds.  Valid for Landsat 8-9 images, and, for Sentinel-2 images with '
-    'the `qa` :option:`--mask-method`.'
+    'the `qa` ``--mask-method``.'
 )
 @click.option(
     '-ms/-nms', '--mask-shadows/--no-mask-shadows', default=True, show_default=True,
@@ -266,7 +291,7 @@ def cli(ctx, verbose, quiet):
 )
 @click.option(
     '-p', '--prob', type=click.FloatRange(min=0, max=100), default=60, show_default=True,
-    help='Cloud probability threshold (%). Valid for Sentinel-2 images with the `cloud-prob` :option:`--mask-method`'
+    help='Cloud probability threshold (%). Valid for Sentinel-2 images with the `cloud-prob` ``--mask-method``'
 )
 @click.option(
     '-d', '--dark', type=click.FloatRange(min=0, max=1), default=.15, show_default=True,
@@ -371,7 +396,7 @@ def config(ctx, mask_cirrus, mask_shadows, mask_method, prob, dark, shadow_dist,
 @click.option(
     '-cp', '--cloudless-portion', type=click.FloatRange(min=0, max=100), default=0, show_default=True,
     help='Lower limit on the cloud/shadow free portion of the region (%).  If cloud/shadow masking is not supported '
-    'for the specified collection, :option:`--cloudless-portion` will operate like :option:`--fill-portion`.'
+    'for the specified collection, ``--cloudless-portion`` will operate like ``--fill-portion``.'
 )
 @click.option(
     '-cf', '--custom-filter', type=click.STRING, default=None,
@@ -474,18 +499,21 @@ def search(
 # download command
 @click.command(cls=ChainedCommand)
 @click.option('-i', '--id', 'image_id', type=click.STRING, multiple=True, help='Earth Engine image ID(s) to download.')
+@crs_option
 @bbox_option
 @region_option
-@click.option(
-    '-dd', '--download-dir', type=click.Path(exists=True, file_okay=False, dir_okay=True, writable=True), default=None,
-    show_default='current working directory.', help='Directory to download image file(s) into.'
-)
-@crs_option
 @scale_option
+@crs_transform_option
+@shape_option
+@like_option
 @dtype_option
 @mask_option
 @resampling_option
 @scale_offset_option
+@click.option(
+    '-dd', '--download-dir', type=click.Path(exists=True, file_okay=False, dir_okay=True, writable=True), default=None,
+    show_default='current working directory.', help='Directory to download image file(s) into.'
+)
 @click.option(
     '-mts', '--max-tile-size', type=click.FLOAT, default=BaseImage._ee_max_tile_size, show_default=True,
     help='Maximum download tile size (MB).'
@@ -496,16 +524,16 @@ def search(
 )
 @click.option('-o', '--overwrite', is_flag=True, default=False, help='Overwrite the destination file if it exists.')
 @click.pass_obj
-def download(obj, image_id, bbox, region, download_dir, mask, max_tile_size, max_tile_dim, overwrite, **kwargs):
+def download(obj, image_id, bbox, region, like, download_dir, mask, max_tile_size, max_tile_dim, overwrite, **kwargs):
     # @formatter:off
     """
     Download image(s).
 
     Download Earth Engine image(s) to GeoTIFF file(s), allowing optional region of interest, and other image
-    formatting options to be specified.  Images larger than the `Earth Engine size limit
-    <https://developers.google.com/earth-engine/apidocs/ee-image-getdownloadurl>`_ are split and downloaded as separate
-    tiles, then re-assembled into a single GeoTIFF.  Downloaded image files are populated with metadata from the Earth
-    Engine image and STAC.
+    formatting options to be specified.  Images larger than the
+    `Earth Engine size limit <https://developers.google.com/earth-engine/apidocs/ee-image-getdownloadurl>`_ are split
+    and downloaded as separate tiles, then re-assembled into a single GeoTIFF.  Downloaded image files are populated
+    with metadata from the Earth Engine image and STAC.
 
     This command can be chained after the ``composite`` command, to download the composite image.  It can also be
     chained after the ``search`` command, in which case the search result images will be downloaded, without the need
@@ -526,8 +554,13 @@ def download(obj, image_id, bbox, region, download_dir, mask, max_tile_size, max
 
     Images from other collections, will contain the FILL_MASK band only.
 
-    If neither ``--bbox`` or ``--region`` are specified, the entire image granule
-    will be downloaded.
+    Bounds and resolution of the downloaded image can be specified with ``--region`` / ``--bbox`` and ``--scale`` /
+    ``--shape``, or ``--crs-transform`` and ``--shape``.  The ``--like`` option will automatically derive ``--crs``,
+    ``--crs-transform`` and ``--shape`` from a provided template raster.  If no bounds are specified (with either
+    ``--region``, or ``--crs-transform`` & ``--shape``), the entire image granule is downloaded.
+
+    When ``--crs``, ``--scale``, ``--crs-transform`` and ``--shape`` are not specified, the pixel grids of the
+    downloaded and Earth Engine images will coincide.
 
     Image filenames are derived from their Earth Engine ID.
     \b
@@ -561,23 +594,26 @@ def download(obj, image_id, bbox, region, download_dir, mask, max_tile_size, max
 # export command
 @click.command(cls=ChainedCommand)
 @click.option('-i', '--id', 'image_id', type=click.STRING, multiple=True, help='Earth Engine image ID(s) to export.')
+@crs_option
 @bbox_option
 @region_option
-@click.option(
-    '-df', '--drive-folder', type=click.STRING, default=None, show_default='root folder.',
-    help='Google Drive folder to export image(s) to.'
-)
-@crs_option
 @scale_option
+@crs_transform_option
+@shape_option
+@like_option
 @dtype_option
 @mask_option
 @resampling_option
 @scale_offset_option
+@click.option(
+    '-df', '--drive-folder', type=click.STRING, default=None, show_default='root folder.',
+    help='Google Drive folder to export image(s) to.'
+)
 @click.option(
     '-w/-nw', '--wait/--no-wait', default=True, show_default=True, help='Whether to wait for the export to complete.'
 )
 @click.pass_obj
-def export(obj, image_id, bbox, region, drive_folder, mask, wait, **kwargs):
+def export(obj, image_id, bbox, region, like, drive_folder, mask, wait, **kwargs):
     # @formatter:off
     """
     Export image(s) to Google Drive.
@@ -604,8 +640,13 @@ def export(obj, image_id, bbox, region, drive_folder, mask, wait, **kwargs):
 
     Images from other collections, will contain the FILL_MASK band only.
 
-    If neither ``--bbox`` or ``--region`` are specified, the entire image granule
-    will be downloaded.
+    Bounds and resolution of the exported image can be specified with ``--region`` / ``--bbox`` and ``--scale`` /
+    ``--shape``, or ``--crs-transform`` and ``--shape``.  The ``--like`` option will automatically derive ``--crs``,
+    ``--crs-transform`` and ``--shape`` from a provided template raster.  If no bounds are specified (with either
+    ``--region``, or ``--crs-transform`` & ``--shape``), the entire image granule is exported.
+
+    When ``--crs``, ``--scale``, ``--crs-transform`` and ``--shape`` are not specified, the pixel grids of the
+    exported and Earth Engine images will coincide.
 
     Image filenames are derived from their Earth Engine ID.
     \b
@@ -662,17 +703,17 @@ def export(obj, image_id, bbox, region, drive_folder, mask, wait, **kwargs):
 @click.option(
     '-b', '--bbox', type=click.FLOAT, nargs=4, default=None, callback=_bbox_cb,
     help='Give preference to images with the highest cloudless (or filled) portion inside this bounding box (left, '
-    'bottom, right, top).  Valid for `mosaic` and `q-mosaic` compositing :option:`--method`.'
+    'bottom, right, top).  Valid for `mosaic` and `q-mosaic` compositing ``--method``.'
 )
 @click.option(
     '-r', '--region', type=click.Path(exists=True, dir_okay=False, allow_dash=True), default=None, callback=_region_cb,
     help='Give preference to images with the highest cloudless (or filled) portion inside this geojson polygon, '
-    'or raster file, region.  Valid for `mosaic` and `q-mosaic` compositing :option:`--method`.'
+    'or raster file, region.  Valid for `mosaic` and `q-mosaic` compositing ``--method``.'
 )
 @click.option(
     '-d', '--date', type=click.DateTime(),
     help='Give preference to images closest to this date (UTC).  Valid for `mosaic` and `q-mosaic` compositing '
-    ':option:`--method`.'
+    '``--method``.'
 )
 @click.pass_obj
 def composite(obj, image_id, mask, method, resampling, bbox, region, date):
diff --git a/geedim/download.py b/geedim/download.py
index 257977c..7d2384b 100644
--- a/geedim/download.py
+++ b/geedim/download.py
@@ -31,7 +31,7 @@
 import rasterio as rio
 from rasterio.crs import CRS
 from rasterio.enums import Resampling as RioResampling
-from rasterio.windows import Window
+from rasterio import windows, features, warp
 from tqdm import TqdmWarning
 from tqdm.auto import tqdm
 from tqdm.contrib.logging import logging_redirect_tqdm
@@ -221,6 +221,14 @@ def refl_bands(self) -> Optional[List[str]]:
             return None
         return [bname for bname, bdict in self._stac.band_props.items() if 'center_wavelength' in bdict]
 
+    @property
+    def profile(self) -> Dict:
+        """ Rasterio image profile. """
+        return dict(
+            crs=self.crs, transform=self.transform, width=self.shape[1], height=self.shape[0], dtype=self.dtype,
+            count=self.count
+        )
+
     @staticmethod
     def _get_projection(ee_info: Dict, min_scale=True) -> Dict:
         """
@@ -364,23 +372,39 @@ def _scale_offset(ee_image: ee.Image, band_properties: List[Dict]) -> ee.Image:
         return ee.Image(adj_im.copyProperties(ee_image, ee_image.propertyNames()))
 
     def _prepare_for_export(
-        self, region: Dict = None, crs: str = None, scale: float = None,
-        resampling: ResamplingMethod = _default_resampling, dtype: str = None, scale_offset: bool = False
+        self, crs: str = None, crs_transform: Tuple[float] = None, shape: Tuple[int, int] = None, region: Dict =
+        None, scale: float = None, resampling: ResamplingMethod = _default_resampling, dtype: str = None,
+        scale_offset: bool = False
     ) -> 'BaseImage':
         """
         Prepare the encapsulated image for export/download.  Will reproject, resample, clip and convert the image
         according to the provided parameters.
 
+        Export bounds and resolution can be specified with ``region`` and ``scale`` / ``shape``, or ``crs_transform``
+        and ``shape``.  If no bounds are specified (with either ``region``, or ``crs_transform`` & ``shape``), the
+        entire image granule is exported.
+
+        When ``crs`` and ``scale`` are not specified (or match those of the encapsulated image), and ``crs_transform``
+        & ``shape`` are not specified, the pixel grid of the exported image will coincide with that of the encapsulated
+        image.
+
         Parameters
         ----------
-        region : dict, geojson, ee.Geometry, optional
-            Region of interest (WGS84) to export.  Defaults tp export the entire image if it has a footprint.
         crs : str, optional
-            WKT, EPSG etc specification of CRS to export to.  Where image bands have different CRSs, all are
+            WKT, EPSG etc. specification of CRS to export to.  Where image bands have different CRSs, all are
             re-projected to this CRS. Defaults to use the CRS of the minimum scale band if available.
+        crs_transform: tuple of float, list of float, rio.Affine, optional
+            List of 6 numbers specifying an affine transform in the specified CRS.  In row-major order:
+            [xScale, xShearing, xTranslation, yShearing, yScale, yTranslation].  All bands are re-projected to
+            this transform.
+        shape: tuple of int, optional
+            (height, width) dimensions to export (pixels).
+        region : dict, geojson, ee.Geometry, optional
+            Region of interest (WGS84) to export.  Defaults to the image footprint, when available.
         scale : float, optional
             Pixel scale (m) to export to.  Where image bands have different scales, all are re-projected to this scale.
-            Defaults to use the minimum scale of image bands if available.
+            Ignored if ``crs`` and ``crs_transform`` are specified.  Defaults to use the minimum scale of
+            image bands if available.
         resampling : ResamplingMethod, optional
             Resampling method - see :class:`~geedim.enums.ResamplingMethod` for available options.
         dtype: str, optional
@@ -395,26 +419,31 @@ def _prepare_for_export(
             Prepared image.
         """
 
-        if not region or not crs or not scale:
-            # One or more of region, crs and scale were not provided, so get the image values to use instead
+        # Check for parameter combination errors.
+        # This necessarily duplicates some of what is done in ee.Image.getDownloadURL, so that errors are raised
+        # before tile creation & download.
+        if (not crs or not region or not scale) and (not crs or not crs_transform or not shape):
             if not self.has_fixed_projection:
-                # Raise an error if this image has no fixed projection
+                # if the image has no fixed projection, either crs, region, & scale; or crs, crs_transform and shape
+                # must be specified
                 raise ValueError(
-                    f'This image does not have a fixed projection, you need to specify a region, '
-                    f'crs and scale.'
+                    f'This image does not have a fixed projection, you need to specify a crs, region & scale; or a '
+                    f'crs, crs_transform & shape.'
                 )
 
-        if not region and not self.footprint:
-            raise ValueError(f'This image does not have a footprint, you need to specify a region.')
-
-        if self.crs == 'EPSG:4326' and not scale:
-            # ee.Image.prepare_for_export() expects a scale in meters, but if the image is EPSG:4326, the default scale
-            # is in degrees.
-            raise ValueError(f'This image is in EPSG:4326, you need to specify a scale in meters.')
+        if (not self.footprint) and (not region and (not crs or not crs_transform or not shape)):
+            # if the image has no footprint (i.e. it is 'unbounded'), either region; or crs, crs_transform and shape
+            # must be specified
+            raise ValueError(
+                f'This image is unbounded, you need to specify a region; or a crs, crs_transform and '
+                f'shape.'
+            )
 
-        region = region or self.footprint
-        crs = crs or ee.Projection(self.crs, tuple(self.transform)[:6])
-        scale = scale or self.scale
+        if self.crs == 'EPSG:4326' and not scale and not shape:
+            # If the image is in EPSG:4326, either scale (in meters); or shape must be specified.
+            # Note that ee.Image.prepare_for_export() expects a scale in meters, but if the image is EPSG:4326,
+            # the default scale is in degrees.
+            raise ValueError(f'This image is in EPSG:4326, you need to specify a scale (in meters); or a shape.')
 
         if crs == 'SR-ORG:6974':
             raise ValueError(
@@ -422,6 +451,11 @@ def _prepare_for_export(
                 'https://issuetracker.google.com/issues/194561313'
             )
 
+        if scale and shape:
+            # This error is raised in later calls to ee.Image.getInfo(), but is neater to raise here first
+            raise ValueError('You can specify one of scale or shape, but not both.')
+
+        # perform image scale/offset, dtype and resampling operations
         ee_image = self.ee_image
         if scale_offset:
             ee_image = self._scale_offset(ee_image, self.band_properties)
@@ -439,10 +473,44 @@ def _prepare_for_export(
             ee_image = utils.resample(ee_image, resampling)
 
         ee_image = self._convert_dtype(ee_image, dtype=dtype or im_dtype)
-        # TODO: Specify `crs_transform` and `dimensions` (as in tile), so that everything stays on the source grid
-        #  where possible i.e. where the export CRS and scale are the same as the source.
-        export_args = dict(region=region, crs=crs, scale=scale, fileFormat='GeoTIFF', filePerBand=False)
-        ee_image, _ = ee_image.prepare_for_export(export_args)
+
+        # configure the export parameter values
+        crs = crs or self.crs
+        if not crs_transform and not shape:
+            # if none of crs_transform, shape, region or scale are specified, then set region and scale to defaults
+            region = region or self.footprint
+            scale = scale or self.scale
+
+            if (crs == self.crs) and (scale == self.scale):
+                # if crs_transform & shape are not already specified, and crs & scale match this image's crs & scale,
+                # then set crs_transform and shape to export on this image's (self) pixel grid
+                if region == self.footprint:
+                    crs_transform = self.transform
+                    shape = self.shape
+                else:
+                    # find a crs_transform and shape that encompasses region
+                    region_crs = region['crs']['properties']['name'] if 'crs' in region else 'EPSG:4326'
+                    region_bounds = warp.transform_bounds(region_crs, crs, *features.bounds(region))
+                    region_win = windows.from_bounds(*region_bounds, transform=self.transform)
+                    region_win = utils.expand_window_to_grid(region_win)
+                    crs_transform = self.transform * rio.Affine.translation(region_win.col_off, region_win.row_off)
+                    shape = (region_win.height, region_win.width)
+
+                # prevent exporting with region & scale now that crs_transform and shape are set
+                region = None
+                scale = None
+
+        # create the export parameter dict
+        crs_transform = crs_transform[:6] if crs_transform else None
+        dimensions = shape[::-1] if shape else None
+        export_kwargs = dict(
+            crs=crs, crs_transform=crs_transform, dimensions=dimensions, region=region, scale=scale,
+            fileFormat='GeoTIFF', filePerBand=False
+        )
+        # drop items with values==None
+        export_kwargs = {k: v for k, v in export_kwargs.items() if v is not None}
+        logger.debug(f'ee.Image.prepare_for_export() params: {export_kwargs}')
+        ee_image, _ = ee_image.prepare_for_export(export_kwargs)
         return BaseImage(ee_image)
 
     def _prepare_for_download(self, set_nodata: bool = True, **kwargs) -> Tuple['BaseImage', Dict]:
@@ -594,7 +662,7 @@ def _tiles(exp_image: 'BaseImage', tile_shape: Tuple[int, int]) -> Iterator[Tile
         for tile_start in start_range:
             tile_stop = np.clip(np.add(tile_start, tile_shape), a_min=None, a_max=image_shape)
             clip_tile_shape = (tile_stop - tile_start).tolist()  # tolist is just to convert to native int
-            tile_window = Window(tile_start[1], tile_start[0], clip_tile_shape[1], clip_tile_shape[0])
+            tile_window = windows.Window(tile_start[1], tile_start[0], clip_tile_shape[1], clip_tile_shape[0])
             yield Tile(exp_image, tile_window)
 
     @staticmethod
@@ -639,6 +707,14 @@ def export(self, filename: str, folder: str = None, wait: bool = True, **kwargs)
         """
         Export the encapsulated image to Google Drive.
 
+        Export bounds and resolution can be specified with ``region`` and ``scale`` / ``shape``, or ``crs_transform``
+        and ``shape``.  If no bounds are specified (with either ``region``, or ``crs_transform`` & ``shape``), the
+        entire image granule is exported.
+
+        When ``crs`` and ``scale`` are not specified (or match those of the encapsulated image), and ``crs_transform``
+        & ``shape`` are not specified, the pixel grid of the exported image will coincide with that of the encapsulated
+        image.
+
         Parameters
         ----------
         filename : str
@@ -647,14 +723,21 @@ def export(self, filename: str, folder: str = None, wait: bool = True, **kwargs)
             Google Drive folder to export to. Defaults to the root.
         wait : bool
             Wait for the export to complete before returning.
-        region : dict, ee.Geometry, optional
-            Region defined by geojson polygon in WGS84. Defaults to the entire image granule.
         crs : str, optional
-            Reproject image(s) to this EPSG or WKT CRS.  Where image bands have different CRSs, all are
-            re-projected to this CRS. Defaults to the CRS of the minimum scale band.
+            WKT, EPSG etc specification of CRS to export to.  Where image bands have different CRSs, all are
+            re-projected to this CRS. Defaults to use the CRS of the minimum scale band if available.
+        crs_transform: tuple of float, list of float, rio.Affine, optional
+            List of 6 numbers specifying an affine transform in the specified CRS.  In row-major order:
+            [xScale, xShearing, xTranslation, yShearing, yScale, yTranslation].  All bands are re-projected to
+            this transform.
+        shape: tuple of int, optional
+            (height, width) dimensions to export (pixels).
+        region : dict, geojson, ee.Geometry, optional
+            Region of interest (WGS84) to export.  Defaults to the image footprint, when available.
         scale : float, optional
-            Resample image(s) to this pixel scale (size) (m).  Where image bands have different scales,
-            all are resampled to this scale. Defaults to the minimum scale of image bands.
+            Pixel scale (m) to export to.  Where image bands have different scales, all are re-projected to this scale.
+            Ignored if ``crs`` and ``crs_transform`` are specified.  Defaults to use the minimum scale of
+            image bands if available.
         resampling : ResamplingMethod, optional
             Resampling method - see :class:`~geedim.enums.ResamplingMethod` for available options.
         dtype: str, optional
@@ -695,6 +778,14 @@ def download(
         separate tiles, then re-assembled into a single GeoTIFF.  Downloaded image files are populated with metadata
         from the Earth Engine image and STAC.
 
+        Image bounds and resolution can be specified with ``region`` and ``scale`` / ``shape``, or ``crs_transform``
+        and ``shape``.  If no bounds are specified (with either ``region``, or ``crs_transform`` & ``shape``), the
+        entire image granule is downloaded.
+
+        When ``crs`` and ``scale`` are not specified (or match those of the encapsulated image), and ``crs_transform``
+        & ``shape`` are not specified, the pixel grid of the downloaded image will coincide with that of the
+        encapsulated image.
+
         Parameters
         ----------
         filename: pathlib.Path, str
@@ -707,16 +798,23 @@ def download(
             Maximum tile size (MB).  If None, defaults to the Earth Engine download size limit (32 MB).
         max_tile_dim: int, optional
             Maximum tile width/height (pixels).  If None, defaults to Earth Engine download limit (10000).
-        region : dict, ee.Geometry, optional
-            Region defined by geojson polygon in WGS84.  Defaults to the entire image granule.
         crs : str, optional
-            Reproject image(s) to this EPSG or WKT CRS.  Where image bands have different CRSs, all are
-            re-projected to this CRS.  Defaults to the CRS of the minimum scale band.
+            WKT, EPSG etc. specification of CRS to export to.  Where image bands have different CRSs, all are
+            re-projected to this CRS. Defaults to use the CRS of the minimum scale band if available.
+        crs_transform: tuple of float, list of float, rio.Affine, optional
+            List of 6 numbers specifying an affine transform in the specified CRS.  In row-major order:
+            [xScale, xShearing, xTranslation, yShearing, yScale, yTranslation].  All bands are re-projected to
+            this transform.
+        shape: tuple of int, optional
+            (height, width) dimensions to export (pixels).
+        region : dict, geojson, ee.Geometry, optional
+            Region of interest (WGS84) to export.  Defaults to the image footprint, when available.
         scale : float, optional
-            Resample image(s) to this pixel scale (size) (m).  Where image bands have different scales,
-            all are resampled to this scale.  Defaults to the minimum scale of image bands.
+            Pixel scale (m) to export to.  Where image bands have different scales, all are re-projected to this scale.
+            Ignored if ``crs`` and ``crs_transform`` are specified.  Defaults to use the minimum scale of
+            image bands if available.
         resampling : ResamplingMethod, optional
-            Resampling method - see :class:`~geedim.enums.ResamplingMethod` for available options.  Defaults to `near`.
+            Resampling method - see :class:`~geedim.enums.ResamplingMethod` for available options.
         dtype: str, optional
            Convert to this data type (`uint8`, `int8`, `uint16`, `int16`, `uint32`, `int32`, `float32` or
            `float64`). Defaults to auto select a minimal type that can represent the range of pixel values.
diff --git a/geedim/utils.py b/geedim/utils.py
index 0471c92..3fe3a8d 100644
--- a/geedim/utils.py
+++ b/geedim/utils.py
@@ -29,9 +29,10 @@
 import requests
 from geedim.enums import ResamplingMethod
 from rasterio.warp import transform_geom
+from rasterio.env import GDALVersion
+from rasterio.windows import Window
 from requests.adapters import HTTPAdapter, Retry
 import numpy as np
-from rasterio.env import GDALVersion
 from tqdm.auto import tqdm
 
 if '__file__' in globals():
@@ -308,3 +309,27 @@ def retry_session(
     return session
 
 
+def expand_window_to_grid(win: Window, expand_pixels: Tuple[int, int] = (0, 0)) -> Window:
+    """
+    Expand rasterio window extents to the nearest whole numbers i.e. for ``expand_pixels`` >= (0, 0), it will return a
+    window that contains the original extents.
+
+    Parameters
+    ----------
+    win: rasterio.windows.Window
+        Window to expand.
+    expand_pixels: tuple, optional
+        Tuple specifying the number of (rows, columns) pixels to expand the window by.
+
+    Returns
+    -------
+    rasterio.windows.Window
+        Expanded window.
+    """
+    col_off, col_frac = np.divmod(win.col_off - expand_pixels[1], 1)
+    row_off, row_frac = np.divmod(win.row_off - expand_pixels[0], 1)
+    width = np.ceil(win.width + 2 * expand_pixels[1] + col_frac)
+    height = np.ceil(win.height + 2 * expand_pixels[0] + row_frac)
+    exp_win = Window(int(col_off), int(row_off), int(width), int(height))
+    return exp_win
+


Test patch:
diff --git a/tests/test_cli.py b/tests/test_cli.py
index f9e148f..d74003a 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -17,7 +17,7 @@
 import pathlib
 from datetime import datetime
 from glob import glob
-from typing import List, Dict
+from typing import List, Dict, Tuple
 
 import ee
 import numpy as np
@@ -30,6 +30,7 @@
 from rasterio.crs import CRS
 from rasterio.features import bounds
 from rasterio.warp import transform_geom
+from rasterio.transform import Affine
 
 
 @pytest.fixture
@@ -89,7 +90,7 @@ def gedi_image_id_list() -> List[str]:
 
 def _test_downloaded_file(
     filename: pathlib.Path, region: Dict = None, crs: str = None, scale: float = None, dtype: str = None,
-    scale_offset: bool = None
+    scale_offset: bool = None, transform: Affine = None, shape: Tuple[int, int] = None
 ):
     """ Helper function to test image file format against given parameters. """
     with rio.open(filename, 'r') as ds:
@@ -119,6 +120,10 @@ def _test_downloaded_file(
             array = ds.read(refl_bands, masked=True)
             assert all(array.min(axis=(1, 2)) >= -0.5)
             assert all(array.max(axis=(1, 2)) <= 1.5)
+        if transform:
+            assert ds.transform[:6] == transform[:6]
+        if shape:
+            assert ds.shape == tuple(shape)
 
 
 @pytest.mark.parametrize(
@@ -304,7 +309,7 @@ def test_search_custom_filter_l9(region_25ha_file: pathlib.Path, runner: CliRunn
         ('gedi_cth_image_id', 'region_25ha_file'),
     ]
 )  # yapf: disable
-def test_download_defaults(
+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.  """
@@ -320,7 +325,68 @@ def test_download_defaults(
     # test downloaded file readability and format
     with open(region_file) as f:
         region = json.load(f)
-    _test_downloaded_file(out_file, region)
+    _test_downloaded_file(out_file, region=region)
+
+
+@pytest.mark.parametrize(
+    '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. """
+    image_id = request.getfixturevalue(image_id)
+    region_file = request.getfixturevalue(region_file)
+    out_file = tmp_path.joinpath(image_id.replace('/', '-') + '.tif')
+
+    # find a transform and shape for region_file
+    with open(region_file) as f:
+        region = json.load(f)
+    region_bounds = bounds(region)
+    crs = 'EPSG:4326'
+    shape = (11, 12)
+    shape_str = ' '.join(map(str, shape))
+    crs_transform = rio.transform.from_bounds(*region_bounds, *shape[::-1])
+    crs_transform_str = ' '.join(map(str, crs_transform[:6]))
+
+    # 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())
+
+    # test downloaded file readability and format
+    _test_downloaded_file(out_file, crs=crs, region=region, transform=crs_transform, shape=shape)
+
+
+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. """
+    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
+    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())
+
+    # 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())
+
+    # 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:
+        assert l8_im.crs == s2_im.crs
+        assert l8_im.shape == s2_im.shape
+        assert l8_im.transform[:6] == s2_im.transform[:6]
 
 
 @pytest.mark.parametrize(
@@ -338,7 +404,7 @@ def test_download_params(
     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 possible cli params. """
+    """ 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')
@@ -419,7 +485,7 @@ def test_composite_defaults(
         ('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),     # TODO revert to GEDI medoid
+        ('gedi_image_id_list', 'medoid', None, None, True, 'bilinear', 25),
     ]
 )
 def test_composite_params(
diff --git a/tests/test_collection.py b/tests/test_collection.py
index 3c2587e..731a0f3 100644
--- a/tests/test_collection.py
+++ b/tests/test_collection.py
@@ -483,7 +483,7 @@ def test_composite_landsat_cloud_mask_params(l8_9_image_list, region_10000ha):
         ('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, {}),    # TODO: include medoid & median?
+        ('gedi_image_list', CompositeMethod.medoid, True, None, None, {}),
     ]
 )  # yapf: disabl;e
 def test_composite(image_list, method, mask, region, date, cloud_kwargs, request):
diff --git a/tests/test_download.py b/tests/test_download.py
index 9cf379e..bb0beff 100644
--- a/tests/test_download.py
+++ b/tests/test_download.py
@@ -23,11 +23,7 @@
 import rasterio as rio
 from geedim.download import BaseImage
 from geedim.enums import ResamplingMethod
-from rasterio import Affine
-from rasterio.coords import BoundingBox
-from rasterio.features import bounds
-from rasterio.warp import transform_geom
-from rasterio.windows import union
+from rasterio import Affine, windows
 
 
 class BaseImageLike:
@@ -70,6 +66,10 @@ def l9_base_image(l9_image_id: str) -> BaseImage:
     """ A BaseImage instance encapsulating a Landsat-9 image.  Covers `region_*ha`.  """
     return BaseImage.from_id(l9_image_id)
 
+@pytest.fixture(scope='session')
+def l8_base_image(l8_image_id: str) -> BaseImage:
+    """ A BaseImage instance encapsulating a Landsat-8 image.  Covers `region_*ha`.  """
+    return BaseImage.from_id(l8_image_id)
 
 @pytest.fixture(scope='session')
 def landsat_ndvi_base_image(landsat_ndvi_image_id: str) -> BaseImage:
@@ -83,6 +83,48 @@ def modis_nbar_base_image(modis_nbar_image_id: str, region_100ha: Dict) -> BaseI
     return BaseImage(ee.Image(modis_nbar_image_id).clip(region_100ha).reproject(crs='EPSG:3857', scale=500))
 
 
+def bounds_polygon(left: float, bottom: float, right: float, top: float, crs: str = None):
+    """ Return a geojson polygon of the given bounds. """
+    coordinates = [[
+        [left, bottom],
+        [left, top],
+        [right, top],
+        [right, bottom],
+        [left, bottom],
+    ]]
+    poly = dict(type='Polygon', coordinates=coordinates)
+    if crs and crs != 'EPSG:4326':
+        poly.update(crs=dict(type='name', properties=dict(name=crs)))
+    return poly
+
+
+def _test_export_profile(exp_profile, tgt_profile, crs_transform=False):
+    """
+    Test an export/download rasterio profile against a target rasterio profile.
+    ``crs_transform`` specifies whether the export/download profile was generated by specifying crs_transform.
+    """
+
+    compare_fn = (lambda x: x) if crs_transform else (lambda x: pytest.approx(x, rel=0.05))
+
+    for key in ['crs', 'dtype', 'count']:
+        assert exp_profile[key] == tgt_profile[key]
+
+    assert exp_profile['transform'][0] == tgt_profile['transform'][0]
+
+    for key in ['width', 'height']:
+        if crs_transform:
+            assert exp_profile[key] == tgt_profile[key]
+    assert exp_profile['transform'][:6] == compare_fn(tgt_profile['transform'][:6])
+
+    tgt_bounds = windows.bounds(
+        windows.Window(0, 0, tgt_profile['width'], tgt_profile['height']), tgt_profile['transform']
+    )
+    exp_bounds = windows.bounds(
+        windows.Window(0, 0, exp_profile['width'], exp_profile['height']), exp_profile['transform']
+    )
+    assert exp_bounds == compare_fn(tgt_bounds)
+
+
 def test_id_name(user_base_image: BaseImage, s2_sr_base_image: BaseImage):
     """ Test `id` and `name` properties for different scenarios. """
     assert user_base_image.id is None
@@ -187,14 +229,52 @@ def test_str_format_size(size: int, exp_str: str):
     assert BaseImage._str_format_size(size) == exp_str
 
 
+@pytest.mark.parametrize(
+    'params', [
+        dict(),
+        dict(crs='EPSG:4326'),
+        dict(crs='EPSG:4326', region='region_25ha'),
+        dict(crs='EPSG:4326', scale=100),
+        dict(crs='EPSG:4326', crs_transform=Affine.identity()),
+        dict(crs='EPSG:4326', shape=(100, 100)),
+        dict(region='region_25ha', scale=100),
+        dict(crs_transform=Affine.identity(), shape=(100, 100)),
+    ]
+)  # yapf: disable
+def test_prepare_no_fixed_projection(user_base_image: BaseImage, params: Dict, request: pytest.FixtureRequest):
+    """
+    Test BaseImage._prepare_for_export raises an exception when the image has no fixed projection, and insufficient
+    crs & region defining parameters are specified.
+    """
+    if 'region' in params:
+        params['region'] = request.getfixturevalue(params['region'])
+    with pytest.raises(ValueError) as ex:
+        user_base_image._prepare_for_export(**params)
+    assert 'does not have a fixed projection' in str(ex)
+
+
+@pytest.mark.parametrize(
+    'params', [
+        dict(),
+        dict(crs='EPSG:4326'),
+        dict(crs='EPSG:4326', scale=100),
+        dict(crs='EPSG:4326', crs_transform=Affine.identity()),
+        dict(crs='EPSG:4326', shape=(100, 100)),
+        dict(crs_transform=Affine.identity(), shape=(100, 100)),
+    ]
+)  # yapf: disable
+def test_prepare_unbounded(user_fix_base_image: BaseImage, params: Dict):
+    """
+    Test BaseImage._prepare_for_export raises an exception when the image is unbounded, and insufficient bounds
+    defining parameters are specified.
+    """
+    with pytest.raises(ValueError) as ex:
+        user_fix_base_image._prepare_for_export(**params)
+    assert 'This image is unbounded' in str(ex)
+
+
 def test_prepare_exceptions(user_base_image: BaseImage, user_fix_base_image: BaseImage, region_25ha: Dict):
-    """ Test BaseImage._prepare_for_export() error cases. """
-    with pytest.raises(ValueError):
-        # no fixed projection and no region / crs / scale
-        user_base_image._prepare_for_export()
-    with pytest.raises(ValueError):
-        # no footprint or region
-        user_fix_base_image._prepare_for_export()
+    """ Test remaining BaseImage._prepare_for_export() error cases. """
     with pytest.raises(ValueError):
         # EPSG:4326 and no scale
         user_fix_base_image._prepare_for_export(region=region_25ha)
@@ -225,55 +305,59 @@ def test_prepare_exceptions(user_base_image: BaseImage, user_fix_base_image: Bas
         ('modis_nbar_base_image', 'modis_nbar_base_image'),
         ('user_base_image', 'modis_nbar_base_image'),
         ('user_fix_base_image', 'modis_nbar_base_image'),
+        # include the below 2 cases to test where export crs & scale match those of the encapsulated image,
+        # and export region is different to the encapsulated footprint
+        ('l9_base_image', 'l8_base_image'),
+        ('l8_base_image', 'l9_base_image'),
     ]
 )  # yapf: disable
 def test_prepare_for_export(src_image: str, tgt_image: str, request: pytest.FixtureRequest):
-    """ Test BaseImage._prepare_for_export() sets properties of export image as expected.  """
+    """ Test BaseImage._prepare_for_export() generates an export image with the target properties. """
     src_image: BaseImage = request.getfixturevalue(src_image)
     tgt_image: BaseImage = request.getfixturevalue(tgt_image)
     if tgt_image == src_image:
-        exp_image = src_image._prepare_for_export()
+        exp_params_list = [dict()] # test default arguments so crs_transform & shape are set to those of src_image
     else:
-        exp_image = src_image._prepare_for_export(
-            crs=tgt_image.crs, scale=tgt_image.scale, region=tgt_image.footprint, dtype=tgt_image.dtype
-        )
+        exp_params_list = [
+            dict(crs=tgt_image.crs, crs_transform=tgt_image.transform, shape=tgt_image.shape, dtype=tgt_image.dtype),
+            dict(crs=tgt_image.crs, scale=tgt_image.scale, region=tgt_image.footprint, dtype=tgt_image.dtype)
+        ]
+    for exp_params in exp_params_list:
+        crs_transform = 'crs_transform' in exp_params or src_image == tgt_image
+        exp_image = src_image._prepare_for_export(**exp_params)
 
-    assert exp_image.crs == tgt_image.crs
-    assert exp_image.scale == tgt_image.scale
-    assert exp_image.count == src_image.count
-    assert exp_image.dtype == tgt_image.dtype
-
-    # Note that exp_image = ee.Image.prepare_for_export(<tgt_image properties>) resamples and can give an exp_image
-    # with different grid and shape compared to tgt_image.  So just test the region here, as that is what is passed
-    # to EE.
-    exp_region = transform_geom(exp_image.footprint['crs']['properties']['name'], 'EPSG:4326', exp_image.footprint)
-    exp_bounds = bounds(exp_region)
-    tgt_bounds = bounds(tgt_image.footprint)
-    assert exp_bounds == pytest.approx(tgt_bounds, rel=.05)
-    # test exp_bounds contain tgt_bounds
-    assert (
-        (exp_bounds[0] <= tgt_bounds[0]) and (exp_bounds[1] <= tgt_bounds[1]) and (exp_bounds[2] >= tgt_bounds[2]) and
-        (exp_bounds[3] >= tgt_bounds[3])
-    )
+        tgt_profile = tgt_image.profile
+        tgt_profile.update(count=src_image.count)
+        _test_export_profile(exp_image.profile, tgt_profile, crs_transform)
+        assert exp_image.scale == tgt_image.scale
 
 
 @pytest.mark.parametrize(
-    'src_image, tgt_image', [('s2_sr_base_image', 's2_sr_base_image'), ('user_base_image', 's2_sr_base_image'), ]
+    'src_image, tgt_image', [
+        ('s2_sr_base_image', 's2_sr_base_image'),
+        ('user_base_image', 's2_sr_base_image'),
+    ]
 )
-def test_prepare_for_download(src_image: str, tgt_image: str, region_25ha: Dict, request: pytest.FixtureRequest):
+def test_prepare_for_download(
+    src_image: str, tgt_image: str, region_25ha: Dict, request: pytest.FixtureRequest
+):
     """ Test BaseImage._prepare_for_download() sets rasterio profile as expected.  """
     src_image: BaseImage = request.getfixturevalue(src_image)
     tgt_image: BaseImage = request.getfixturevalue(tgt_image)
     if tgt_image == src_image:
-        exp_image, exp_profile = src_image._prepare_for_download()
+        exp_param_list = [dict()]
     else:
-        exp_image, exp_profile = src_image._prepare_for_download(
-            crs=tgt_image.crs, scale=tgt_image.scale, region=region_25ha, dtype=tgt_image.dtype
-        )
-    assert f'EPSG:{exp_profile["crs"].to_epsg()}' == tgt_image.crs
-    assert exp_profile['count'] == src_image.count
-    assert exp_profile['dtype'] == tgt_image.dtype
-    assert exp_profile['nodata'] is not None
+        exp_param_list = [
+            dict(crs=tgt_image.crs, crs_transform=tgt_image.transform, shape=tgt_image.shape, dtype=tgt_image.dtype),
+            dict(crs=tgt_image.crs, scale=tgt_image.scale, region=tgt_image.footprint, dtype=tgt_image.dtype),
+        ]
+
+    for exp_params in exp_param_list:
+        exp_image, exp_profile = src_image._prepare_for_download(**exp_params)
+        tgt_profile = tgt_image.profile
+        tgt_profile.update(count=src_image.count)
+        _test_export_profile(exp_profile, tgt_profile, ('crs_transform' in exp_params) or (tgt_image == src_image))
+        assert exp_profile['nodata'] is not None
 
 
 @pytest.mark.parametrize(
@@ -363,7 +447,7 @@ def test_tiles(image_shape: Tuple, tile_shape: Tuple, image_transform: Affine):
     prev_tile = tiles[0]
     accum_window = prev_tile.window
     for tile in tiles[1:]:
-        accum_window = union(accum_window, tile.window)
+        accum_window = windows.union(accum_window, tile.window)
         assert all(np.array(tile._shape) <= np.array(tile_shape))
 
         if tile.window.row_off == prev_tile.window.row_off:
@@ -380,40 +464,26 @@ def test_tiles(image_shape: Tuple, tile_shape: Tuple, image_transform: Affine):
     assert (accum_window.height, accum_window.width) == exp_image.shape
 
 
-@pytest.mark.parametrize(
-    'base_image, region', [
-        ('user_base_image', 'region_25ha'),
-        ('user_fix_base_image', 'region_25ha'),
-        # ('s2_sr_base_image', 'region_25ha'),
-        # ('l9_base_image', 'region_25ha'),
-        # ('modis_nbar_base_image', 'region_25ha'),
+def test_download(user_base_image: str, tmp_path: pathlib.Path, request: pytest.FixtureRequest):
+    """ Test download file properties and pixel data. """
+    tgt_prof = dict(crs='EPSG:3857', transform=Affine(1, 0, 0, 0, -1, 0), width=10, height=10, dtype='uint16', count=3)
+
+    # form download parameters from tgt_prof
+    shape = (tgt_prof['height'], tgt_prof['width'])
+    tgt_bounds = windows.bounds(
+        windows.Window(0, 0, *shape[::-1]), tgt_prof['transform']
+    )
+    region = bounds_polygon(*tgt_bounds, crs=tgt_prof['crs'])
+    download_params_list = [
+        dict(crs=tgt_prof['crs'], region=region, scale=tgt_prof['transform'][0], dtype=tgt_prof['dtype']),
+        dict(crs=tgt_prof['crs'], crs_transform=tgt_prof['transform'], shape=shape, dtype=tgt_prof['dtype'])
     ]
-)  # yapf: disable
-def test_download(base_image: str, region: Dict, tmp_path: pathlib.Path, request: pytest.FixtureRequest):
-    """ Test downloaded file properties and pixel data.  """
-    base_image = request.getfixturevalue(base_image)
-    region = request.getfixturevalue(region)
-    filename = tmp_path.joinpath('test_user_download.tif')
-    crs = 'EPSG:3857'
-    dtype = 'uint16'
-    scale = 30
-    base_image.download(filename, region=region, crs=crs, scale=scale, dtype=dtype)
-    exp_region = transform_geom('EPSG:4326', crs, region)
-    exp_bounds = BoundingBox(*bounds(exp_region))
-    assert filename.exists()
-    with rio.open(filename, 'r') as ds:
-        assert ds.count == base_image.count
-        assert ds.dtypes[0] == dtype
-        assert ds.nodata == 0
-        assert abs(ds.transform[0]) == scale
-        assert ds.transform.xoff <= exp_bounds.left
-        assert ds.transform.yoff >= exp_bounds.top  # 'EPSG:3857' has y -ve
-        ds_bounds = ds.bounds
-        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])
-        )
-        if ds.count < 4:
+    filename_list = [tmp_path.joinpath('test1.tif'), tmp_path.joinpath('test2.tif')]
+    for filename, download_params in zip(filename_list, download_params_list):
+        user_base_image.download(filename, **download_params)
+        assert filename.exists()
+        with rio.open(filename, 'r') as ds:
+            _test_export_profile(ds.profile, tgt_prof, 'crs_transform' in download_params)
             array = ds.read()
             for i in range(ds.count):
                 assert np.all(array[i] == i + 1)