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

Issue:
BadZipFile: File is not a zip file
I used a library called geemap. One of the functions in this library, `geemap.download_ee_image`, uses the geedim library to download image results calculated through the Earth Engine API.

However, when I use this function to download the image, I can see the progress bar increasing, but when it reaches a certain download progress, it throws an error.



BadZipFile Traceback (most recent call last)
Cell In[23], line 37
34 temp_image.select(["SR_B1", "SR_B2", "SR_B3", "SR_B4","SR_B5", "SR_B7"])
36 y1=1990+i
---> 37 geemap.download_ee_image(temp_image, filename=f'landsat_image_{y1}.tif', scale=30, region=boundary.geometry(),crs='EPSG:4326')
40 for i in range(0,11):
41 start_time=ee.Date('2013-06-01').advance(i,'year')

File d:\anaconda3\envs\myenv\lib\site-packages\geemap\common.py:12925, in download_ee_image(image, filename, region, crs, crs_transform, scale, resampling, dtype, overwrite, num_threads, max_tile_size, max_tile_dim, shape, scale_offset, unmask_value, **kwargs)
12922 kwargs["scale_offset"] = scale_offset
12924 img = gd.download.BaseImage(image)

12925 img.download(filename, overwrite=overwrite, num_threads=num_threads, **kwargs)

File d:\anaconda3\envs\myenv\lib\site-packages\geedim\download.py:955, in BaseImage.download(self, filename, overwrite, num_threads, max_tile_size, max_tile_dim, **kwargs)
953 logger.info(f'Exception: {str(ex)}\nCancelling...')
954 executor.shutdown(wait=False)
--> 955 raise ex
957 bar.update(bar.total - bar.n) # ensure the bar reaches 100%
958 # populate GeoTIFF metadata

File d:\anaconda3\envs\myenv\lib\site-packages\geedim\download.py:951, in BaseImage.download(self, filename, overwrite, num_threads, max_tile_size, max_tile_dim, **kwargs)
949 try:
950 for future in as_completed(futures):
--> 951 future.result()
952 except Exception as ex:
953 logger.info(f'Exception: {str(ex)}\nCancelling...')

File d:\anaconda3\envs\myenv\lib\concurrent\futures_base.py:437, in Future.result(self, timeout)
435 raise CancelledError()
436 elif self._state == FINISHED:
--> 437 return self.__get_result()
439 self._condition.wait(timeout)
441 if self._state in [CANCELLED, CANCELLED_AND_NOTIFIED]:

File d:\anaconda3\envs\myenv\lib\concurrent\futures_base.py:389, in Future.__get_result(self)
387 if self._exception:
388 try:
--> 389 raise self._exception
390 finally:
391 # Break a reference cycle with the exception in self._exception
392 self = None

File d:\anaconda3\envs\myenv\lib\concurrent\futures\thread.py:57, in _WorkItem.run(self)
54 return
56 try:
---> 57 result = self.fn(*self.args, **self.kwargs)
58 except BaseException as exc:
59 self.future.set_exception(exc)

File d:\anaconda3\envs\myenv\lib\site-packages\geedim\download.py:941, in BaseImage.download.<locals>.download_tile(tile)
939 def download_tile(tile):
940 """Download a tile and write into the destination GeoTIFF. """
--> 941 tile_array = tile.download(session=session, bar=bar)
942 with out_lock:
943 out_ds.write(tile_array, window=tile.window)

File d:\anaconda3\envs\myenv\lib\site-packages\geedim\tile.py:121, in Tile.download(self, session, response, bar)
118 zip_buffer.flush()
120 # extract geotiff from zipped buffer into another buffer
--> 121 zip_file = zipfile.ZipFile(zip_buffer)
122 ext_buffer = BytesIO(zip_file.read(zip_file.filelist[0]))
124 # read the geotiff with a rasterio memory file

File d:\anaconda3\envs\myenv\lib\zipfile.py:1269, in ZipFile.init(self, file, mode, compression, allowZip64, compresslevel, strict_timestamps)
1267 try:
1268 if mode == 'r':
-> 1269 self._RealGetContents()
1270 elif mode in ('w', 'x'):
1271 # set the modified flag so central directory gets written
1272 # even if no files are added to the archive
1273 self._didModify = True

File d:\anaconda3\envs\myenv\lib\zipfile.py:1336, in ZipFile._RealGetContents(self)
1334 raise BadZipFile("File is not a zip file")
1335 if not endrec:
-> 1336 raise BadZipFile("File is not a zip file")
1337 if self.debug > 1:
1338 print(endrec)

BadZipFile: File is not a zip file



How can I solve this problem? Is there a better way to directly download Earth Engine image objects?








Discussion / hints:
I briefly reviewed the source code and found that the image download uses a thread pool. However, if one thread encounters an issue, it terminates all download tasks. I think this approach is inefficient in terms of time and data transfer. Is there a way to optimize it? Perhaps adding a retry mechanism to threads that encounter issues could help.
Hello.  The thread pool uses a requests session that is configured for retries.  It is intended to make the download time efficient.

I haven't seen this before.  Can you provide a short code snippet that reproduces the issue please?
I have seen this error before, but have not been able to come up with a consistently reproducible example myself. Often re-trying the download of same extent works as expected.

Could possibly add additional error handling to [tile.py L121](https://github.com/leftfield-geospatial/geedim/blob/fa46cb453a69eac74e7145f44c2e4bf09d688359/geedim/tile.py#L121) in the event that the ZIP file is not correctly downloaded/formatted to re-try the tile? 
Thanks for the info @brownag.  Theoretically, Requests should be doing retries.  I'm not sure why it doesn't raise an error in this situtation.  I need a way of replicating and testing this before implementing anything. 
```python
ee.Initialize()



def rm_cloud(image):
    cloud_shadow_bit_mask = (1 << 3)
    clouds_bit_mask = (1 << 5)
    qa = image.select("QA_PIXEL")
    mask = qa.bitwiseAnd(cloud_shadow_bit_mask).eq(0).And(qa.bitwiseAnd(clouds_bit_mask).eq(0))
    return image.updateMask(mask)


def apply_scale_factors5(image):
    optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
    thermal_bands = image.select('ST_B6').multiply(0.00341802).add(149.0)
    image = image.addBands(optical_bands, None, True)
    image = image.addBands(thermal_bands, None, True)
    return image



def apply_scale_factors8(image):
    optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
    thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
    image = image.addBands(optical_bands, None, True)
    image = image.addBands(thermal_bands, None, True)
    return image



shp_path = r"D:\code_space\gee_test\yyh.shp"
boundary = ee.FeatureCollection(shapefile.Reader(shp_path).__geo_interface__)


if __name__ == '__main__':
    start_time = ee.Date('1990-06-01').advance(1,'year')
    temp_image= ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').filterBounds(boundary.geometry()).filterDate(ee.Date(start_time),ee.Date(start_time).advance((3),'month').advance((29),'day')).map(rm_cloud).map(apply_scale_factors5).median().clip(boundary)
    temp_image.select(["SR_B1", "SR_B2", "SR_B3", "SR_B4","SR_B5", "SR_B7"])  
    y1=1990+1
    geemap.download_ee_image(temp_image, filename=f'./landsat_image_{y1}.tif',  scale=30, region=boundary.geometry(), crs='EPSG:4326')
```



<img width="1223" alt="3bb767ea89f8579d37de753dc4b7700" src="https://github.com/leftfield-geospatial/geedim/assets/59814105/6cad99df-b806-4593-8d19-6f92d9bdf865">




"I don't know if it's because the file I downloaded is too large, or for some other reason?"










```
BadZipFile                                Traceback (most recent call last)
Cell In[1], [line 49](vscode-notebook-cell:?execution_count=1&line=49)
     [47](vscode-notebook-cell:?execution_count=1&line=47) temp_image.select(["SR_B1", "SR_B2", "SR_B3", "SR_B4","SR_B5", "SR_B7"])  
     [48](vscode-notebook-cell:?execution_count=1&line=48) y1=1990+1
---> [49](vscode-notebook-cell:?execution_count=1&line=49) geemap.download_ee_image(temp_image, filename=f'./landsat_image_{y1}.tif',  scale=30, region=boundary.geometry(), crs='EPSG:4326')

File d:\anaconda3\envs\myenv\lib\site-packages\geemap\common.py:12925, in download_ee_image(image, filename, region, crs, crs_transform, scale, resampling, dtype, overwrite, num_threads, max_tile_size, max_tile_dim, shape, scale_offset, unmask_value, **kwargs)
  [12922](file:///D:/anaconda3/envs/myenv/lib/site-packages/geemap/common.py:12922)     kwargs["scale_offset"] = scale_offset
  [12924](file:///D:/anaconda3/envs/myenv/lib/site-packages/geemap/common.py:12924) img = gd.download.BaseImage(image)
> [12925](file:///D:/anaconda3/envs/myenv/lib/site-packages/geemap/common.py:12925) img.download(filename, overwrite=overwrite, num_threads=num_threads, **kwargs)

File d:\anaconda3\envs\myenv\lib\site-packages\geedim\download.py:955, in BaseImage.download(self, filename, overwrite, num_threads, max_tile_size, max_tile_dim, **kwargs)
    [953](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:953)         logger.info(f'Exception: {str(ex)}\nCancelling...')
    [954](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:954)         executor.shutdown(wait=False)
--> [955](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:955)         raise ex
    [957](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:957) bar.update(bar.total - bar.n)   # ensure the bar reaches 100%
    [958](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:958) # populate GeoTIFF metadata

File d:\anaconda3\envs\myenv\lib\site-packages\geedim\download.py:951, in BaseImage.download(self, filename, overwrite, num_threads, max_tile_size, max_tile_dim, **kwargs)
    [949](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:949) try:
    [950](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:950)     for future in as_completed(futures):
--> [951](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:951)         future.result()
    [952](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:952) except Exception as ex:
    [953](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:953)     logger.info(f'Exception: {str(ex)}\nCancelling...')

File d:\anaconda3\envs\myenv\lib\concurrent\futures\_base.py:437, in Future.result(self, timeout)
    [435](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:435)     raise CancelledError()
    [436](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:436) elif self._state == FINISHED:
--> [437](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:437)     return self.__get_result()
    [439](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:439) self._condition.wait(timeout)
    [441](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:441) if self._state in [CANCELLED, CANCELLED_AND_NOTIFIED]:

File d:\anaconda3\envs\myenv\lib\concurrent\futures\_base.py:389, in Future.__get_result(self)
    [387](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:387) if self._exception:
    [388](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:388)     try:
--> [389](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:389)         raise self._exception
    [390](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:390)     finally:
    [391](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:391)         # Break a reference cycle with the exception in self._exception
    [392](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/_base.py:392)         self = None

File d:\anaconda3\envs\myenv\lib\concurrent\futures\thread.py:57, in _WorkItem.run(self)
     [54](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/thread.py:54)     return
     [56](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/thread.py:56) try:
---> [57](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/thread.py:57)     result = self.fn(*self.args, **self.kwargs)
     [58](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/thread.py:58) except BaseException as exc:
     [59](file:///D:/anaconda3/envs/myenv/lib/concurrent/futures/thread.py:59)     self.future.set_exception(exc)

File d:\anaconda3\envs\myenv\lib\site-packages\geedim\download.py:941, in BaseImage.download.<locals>.download_tile(tile)
    [939](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:939) def download_tile(tile):
    [940](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:940)     """Download a tile and write into the destination GeoTIFF. """
--> [941](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:941)     tile_array = tile.download(session=session, bar=bar)
    [942](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:942)     with out_lock:
    [943](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/download.py:943)         out_ds.write(tile_array, window=tile.window)

File d:\anaconda3\envs\myenv\lib\site-packages\geedim\tile.py:121, in Tile.download(self, session, response, bar)
    [118](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/tile.py:118) zip_buffer.flush()
    [120](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/tile.py:120) # extract geotiff from zipped buffer into another buffer
--> [121](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/tile.py:121) zip_file = zipfile.ZipFile(zip_buffer)
    [122](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/tile.py:122) ext_buffer = BytesIO(zip_file.read(zip_file.filelist[0]))
    [124](file:///D:/anaconda3/envs/myenv/lib/site-packages/geedim/tile.py:124) # read the geotiff with a rasterio memory file

File d:\anaconda3\envs\myenv\lib\zipfile.py:1269, in ZipFile.__init__(self, file, mode, compression, allowZip64, compresslevel, strict_timestamps)
   [1267](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1267) try:
   [1268](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1268)     if mode == 'r':
-> [1269](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1269)         self._RealGetContents()
   [1270](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1270)     elif mode in ('w', 'x'):
   [1271](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1271)         # set the modified flag so central directory gets written
   [1272](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1272)         # even if no files are added to the archive
   [1273](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1273)         self._didModify = True

File d:\anaconda3\envs\myenv\lib\zipfile.py:1336, in ZipFile._RealGetContents(self)
   [1334](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1334)     raise BadZipFile("File is not a zip file")
   [1335](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1335) if not endrec:
-> [1336](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1336)     raise BadZipFile("File is not a zip file")
   [1337](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1337) if self.debug > 1:
   [1338](file:///D:/anaconda3/envs/myenv/lib/zipfile.py:1338)     print(endrec)

BadZipFile: File is not a zip file
```
> Thanks for the info @brownag. Theoretically, Requests should be doing retries. I'm not sure why it doesn't raise an error in this situtation. I need a way of replicating and testing this before implementing anything.

Well, I think an unhandled error (i.e. `BadZipFile("File is not a zip file")`) is raised by `zipfile`. The download appears to complete successfully, but the payload is not a valid ZIP file, causing subsequent unpacking of the result to fail. I imagine some sort of atypical response is occurring but not detected. It may be network specific--I have encountered it probably most often when on my work network

Understood that you need a reprex. Will try to find a case where this can happen fairly consistently. Like @mht2953658596 says I have seen this with fairly large downloads--where just one tile throwing this error stops the whole process. I imagine it is not specific to the total download size but rather just that the number of tiles involved makes it more likely that one will fail.

![1abf7b198ebfcd37fd85f98d76bb81d](https://github.com/leftfield-geospatial/geedim/assets/59814105/a6576b6f-7f70-4afc-8b13-a63b779fe23a)

> ```python
> ee.Initialize()
> 
> 
> 
> def rm_cloud(image):
>     cloud_shadow_bit_mask = (1 << 3)
>     clouds_bit_mask = (1 << 5)
>     qa = image.select("QA_PIXEL")
>     mask = qa.bitwiseAnd(cloud_shadow_bit_mask).eq(0).And(qa.bitwiseAnd(clouds_bit_mask).eq(0))
>     return image.updateMask(mask)
> 
> 
> def apply_scale_factors5(image):
>     optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
>     thermal_bands = image.select('ST_B6').multiply(0.00341802).add(149.0)
>     image = image.addBands(optical_bands, None, True)
>     image = image.addBands(thermal_bands, None, True)
>     return image
> 
> 
> 
> def apply_scale_factors8(image):
>     optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
>     thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
>     image = image.addBands(optical_bands, None, True)
>     image = image.addBands(thermal_bands, None, True)
>     return image
> 
> 
> 
> shp_path = r"D:\code_space\gee_test\yyh.shp"
> boundary = ee.FeatureCollection(shapefile.Reader(shp_path).__geo_interface__)
> 
> 
> if __name__ == '__main__':
>     start_time = ee.Date('1990-06-01').advance(1,'year')
>     temp_image= ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').filterBounds(boundary.geometry()).filterDate(ee.Date(start_time),ee.Date(start_time).advance((3),'month').advance((29),'day')).map(rm_cloud).map(apply_scale_factors5).median().clip(boundary)
>     temp_image.select(["SR_B1", "SR_B2", "SR_B3", "SR_B4","SR_B5", "SR_B7"])  
>     y1=1990+1
>     geemap.download_ee_image(temp_image, filename=f'./landsat_image_{y1}.tif',  scale=30, region=boundary.geometry(), crs='EPSG:4326')
> ```
> 
> <img alt="3bb767ea89f8579d37de753dc4b7700" width="1223" src="https://private-user-images.githubusercontent.com/59814105/340699223-6cad99df-b806-4593-8d19-6f92d9bdf865.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MTg3MjcxMzIsIm5iZiI6MTcxODcyNjgzMiwicGF0aCI6Ii81OTgxNDEwNS8zNDA2OTkyMjMtNmNhZDk5ZGYtYjgwNi00NTkzLThkMTktNmY5MmQ5YmRmODY1LnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA2MTglMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwNjE4VDE2MDcxMlomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTVmMjllZmI2MjdjM2VmZmZlYjY5NGRkYWMwNmJmNzc2NzExMmU3OTE4MjJhMTc1MDU1NWY1Y2M4ODRiZTFiMmEmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.wbIF24kpB7o5NWpiKq67CdgFmfMP-4Bjq0TcFL-RY0s">
> 
> "I don't know if it's because the file I downloaded is too large, or for some other reason?"

Thanks for the code @mht2953658596 .  Unfortunately, I can't replicate  - I downloaded it 10x without any issue.  I used the region:

``` python
region = {
    'type': 'Polygon', 'coordinates': [[(9.5, 49), (9.5, 50), (8, 50), (8, 49), (9.5, 49)]]
}
```
as I don't have access to your shape file.  It generates a ~3GB image, similar to yours.

> ![1abf7b198ebfcd37fd85f98d76bb81d](https://private-user-images.githubusercontent.com/59814105/340985508-a6576b6f-7f70-4afc-8b13-a63b779fe23a.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MTg3ODY2MTYsIm5iZiI6MTcxODc4NjMxNiwicGF0aCI6Ii81OTgxNDEwNS8zNDA5ODU1MDgtYTY1NzZiNmYtN2Y3MC00YWZjLThiMTMtYTYzYjc3OWZlMjNhLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA2MTklMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwNjE5VDA4MzgzNlomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPWM3ZDViODYwYWRhOGVmMWZkZDZhNmY4MDU1OTFjYTYzNTZlY2NmMzg0ZWNlYmY2ZmY3ZTQ5NThiNGI5OGRhOGEmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.uNpoaiW6wWHrv0MCN1-JpbneJWSusxDx4LAckp02eWg)

This is more what I was expecting where Requests raises an error.  It does appear to be a network specific reliability issue.  I am going to work on a retry for the tile download that will hopefully fix or at least reduce these problems.

@dugalh  Thank you for your effort. I hope you can add a feature to retry the download of individual tiles in case of timeout or failure. This way, a single tile download failure or blockage won't prevent the entire image from being downloaded.



Fix patch:
diff --git a/geedim/cli.py b/geedim/cli.py
index 4a8fa7a..da1142c 100644
--- a/geedim/cli.py
+++ b/geedim/cli.py
@@ -13,6 +13,7 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
+
 import json
 import logging
 import os
@@ -26,13 +27,14 @@
 import rasterio as rio
 import rasterio.crs as rio_crs
 from click.core import ParameterSource
+from rasterio.errors import CRSError
+
 from geedim import schema, Initialize, version
 from geedim.collection import MaskedCollection
 from geedim.download import BaseImage, supported_dtypes
 from geedim.enums import CloudMaskMethod, CompositeMethod, ResamplingMethod, ExportType
 from geedim.mask import MaskedImage
 from geedim.utils import get_bounds, Spinner, asset_id
-from rasterio.errors import CRSError
 
 logger = logging.getLogger(__name__)
 
@@ -54,18 +56,18 @@ class ChainedCommand(click.Command):
     """
 
     def get_help(self, ctx):
-        """ Strip some RST markup from the help text for CLI display.  Assumes no grid tables. """
+        """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
 
         sub_strings = {
-            '\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'
-            '::': ':',                      # convert from RST '::' to ':'
-            '``(.*?)``': r'\g<1>',          # convert from RST '``literal``' to 'literal'
+            '\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'
+            '::': ':',  # 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'
+            ':option:`(.*?)`': r'\g<1>',  # convert ':option:`--name`' to '--name'
             '`([^<]*) <([^>]*)>`_': r'\g<1>',  # convert from RST cross-ref '`<name> <<link>>`_' to 'name'
         }  # yapf: disable
 
@@ -187,7 +189,7 @@ def _export_type_cb(ctx, param, value):
     return ExportType(value)
 
 
-def _prepare_image_list(obj: SimpleNamespace, mask=False) -> List[MaskedImage, ]:
+def _prepare_image_list(obj: SimpleNamespace, mask=False) -> List[MaskedImage,]:
     """Validate and prepare the obj.image_list for export/download.  Returns a list of MaskedImage objects."""
     if len(obj.image_list) == 0:
         raise click.BadOptionUsage(
@@ -209,54 +211,97 @@ def _prepare_image_list(obj: SimpleNamespace, mask=False) -> List[MaskedImage, ]
 
 # 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,
+    '-b',
+    '--bbox',
+    type=click.FLOAT,
+    nargs=4,
+    default=None,
+    callback=_bbox_cb,
     metavar='LEFT BOTTOM RIGHT TOP',
-    help='Region defined by WGS84 bounding box co-ordinates.'
+    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,
-    help='Region defined by geojson polygon file, or raster file.  Use "-" to read geojson from stdin.'
+    '-r',
+    '--region',
+    type=click.Path(exists=True, dir_okay=False, allow_dash=True),
+    default=None,
+    callback=_region_cb,
+    help='Region defined by geojson polygon file, or raster file.  Use "-" to read geojson from stdin.',
 )
 crs_option = click.option(
-    '-c', '--crs', type=click.STRING, default=None, callback=_crs_cb, show_default='source image CRS.',
-    help='CRS to reproject image(s) to (EPSG string or path to WKT text file).'
+    '-c',
+    '--crs',
+    type=click.STRING,
+    default=None,
+    callback=_crs_cb,
+    show_default='source image CRS.',
+    help='CRS to reproject image(s) to (EPSG string or path to WKT text file).',
 )
 scale_option = click.option(
-    '-s', '--scale', type=click.FLOAT, default=None, show_default='minimum scale of the source image bands.',
-    help='Pixel scale (size) to resample image(s) to (m).'
+    '-s',
+    '--scale',
+    type=click.FLOAT,
+    default=None,
+    show_default='minimum scale of the source image bands.',
+    help='Pixel scale (size) to resample image(s) to (m).',
 )
 dtype_option = click.option(
-    '-dt', '--dtype', type=click.Choice(supported_dtypes, case_sensitive=False), default=None,
+    '-dt',
+    '--dtype',
+    type=click.Choice(supported_dtypes, case_sensitive=False),
+    default=None,
     show_default='smallest data type able to represent the range of pixel values.',
-    help='Data type to convert image(s) to.'
+    help='Data type to convert image(s) to.',
 )
 mask_option = click.option(
-    '-m/-nm', '--mask/--no-mask', default=MaskedImage._default_mask, show_default=True,
+    '-m/-nm',
+    '--mask/--no-mask',
+    default=MaskedImage._default_mask,
+    show_default=True,
     help='Whether to apply cloud/shadow mask(s); or fill mask(s), in the case of images without '
-    'support for cloud/shadow masking.'
+    'support for cloud/shadow masking.',
 )
 resampling_option = click.option(
-    '-rs', '--resampling', type=click.Choice([rm.value for rm in ResamplingMethod], case_sensitive=True),
-    default=BaseImage._default_resampling.value, show_default=True, callback=_resampling_method_cb,
-    help='Resampling method.'
+    '-rs',
+    '--resampling',
+    type=click.Choice([rm.value for rm in ResamplingMethod], case_sensitive=True),
+    default=BaseImage._default_resampling.value,
+    show_default=True,
+    callback=_resampling_method_cb,
+    help='Resampling method.',
 )
 scale_offset_option = click.option(
-    '-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.'
+    '-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,
+    '-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.'
+    '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).'
+    '-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``.'
+    '-l',
+    '--like',
+    type=click.Path(exists=True, dir_okay=False),
+    default=None,
+    help='Template raster from which to derive ``--crs``, ``--crs-transform`` & ``--shape``.',
 )
 
 
@@ -267,7 +312,7 @@ def _prepare_image_list(obj: SimpleNamespace, mask=False) -> List[MaskedImage, ]
 @click.version_option(version=version.__version__, message='%(version)s')
 @click.pass_context
 def cli(ctx, verbose, quiet):
-    """ Search, composite and download Google Earth Engine imagery. """
+    """Search, composite and download Google Earth Engine imagery."""
     ctx.obj = SimpleNamespace(image_list=[], region=None, cloud_kwargs={})
     verbosity = verbose - quiet
     _configure_logging(verbosity)
@@ -279,48 +324,82 @@ def cli(ctx, verbose, quiet):
 # 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'))
+@cli.command(cls=ChainedCommand, context_settings=dict(auto_envvar_prefix='GEEDIM'))
 @click.option(
-    '-mc/-nmc', '--mask-cirrus/--no-mask-cirrus', default=True, show_default=True,
+    '-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` ``--mask-method``.'
+    'the `qa` ``--mask-method``.',
 )
 @click.option(
-    '-ms/-nms', '--mask-shadows/--no-mask-shadows', default=True, show_default=True,
-    help='Whether to mask cloud shadows.'
+    '-ms/-nms',
+    '--mask-shadows/--no-mask-shadows',
+    default=True,
+    show_default=True,
+    help='Whether to mask cloud shadows.',
 )
 @click.option(
-    '-mm', '--mask-method', type=click.Choice([cmm.value for cmm in CloudMaskMethod], case_sensitive=True),
-    default=CloudMaskMethod.cloud_prob.value, show_default=True, callback=_mask_method_cb,
-    help='Method used to mask clouds.  Valid for Sentinel-2 images. '
+    '-mm',
+    '--mask-method',
+    type=click.Choice([cmm.value for cmm in CloudMaskMethod], case_sensitive=True),
+    default=CloudMaskMethod.cloud_prob.value,
+    show_default=True,
+    callback=_mask_method_cb,
+    help='Method used to mask clouds.  Valid for Sentinel-2 images. ',
 )
 @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` ``--mask-method``'
+    '-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` ``--mask-method``',
 )
 @click.option(
-    '-d', '--dark', type=click.FloatRange(min=0, max=1), default=.15, show_default=True,
+    '-d',
+    '--dark',
+    type=click.FloatRange(min=0, max=1),
+    default=0.15,
+    show_default=True,
     help='NIR reflectance threshold for shadow masking. NIR values below this threshold are '
-    'potential cloud shadows.  Valid for Sentinel-2 images'
+    'potential cloud shadows.  Valid for Sentinel-2 images',
 )
 @click.option(
-    '-sd', '--shadow-dist', type=click.INT, default=1000, show_default=True,
-    help='Maximum distance (m) to look for cloud shadows from cloud edges.  Valid for Sentinel-2 images.'
+    '-sd',
+    '--shadow-dist',
+    type=click.INT,
+    default=1000,
+    show_default=True,
+    help='Maximum distance (m) to look for cloud shadows from cloud edges.  Valid for Sentinel-2 images.',
 )
 @click.option(
-    '-b', '--buffer', type=click.INT, default=50, show_default=True,
-    help='Distance (m) to dilate cloud/shadow.  Valid for Sentinel-2 images.'
+    '-b',
+    '--buffer',
+    type=click.INT,
+    default=50,
+    show_default=True,
+    help='Distance (m) to dilate cloud/shadow.  Valid for Sentinel-2 images.',
 )
 @click.option(
-    '-cdi', '--cdi-thresh', type=click.FloatRange(min=-1, max=1), default=None,
+    '-cdi',
+    '--cdi-thresh',
+    type=click.FloatRange(min=-1, max=1),
+    default=None,
     help='Cloud Displacement Index (CDI) threshold.  Values below this threshold are considered potential clouds.  '
-    'Valid for Sentinel-2 images.  By default, the CDI is not used.'
+    'Valid for Sentinel-2 images.  By default, the CDI is not used.',
 )
 @click.option(
-    '-mcd', '--max-cloud-dist', type=click.INT, default=5000, show_default=True,
+    '-mcd',
+    '--max-cloud-dist',
+    type=click.INT,
+    default=5000,
+    show_default=True,
     help='Maximum distance (m) to look for clouds.  Used to form the cloud distance band for the `q-mosaic` '
-    'compositing ``--method``.'
+    'compositing ``--method``.',
 )
 @click.pass_context
 def config(ctx, mask_cirrus, mask_shadows, mask_method, prob, dark, shadow_dist, buffer, cdi_thresh, max_cloud_dist):
@@ -376,55 +455,97 @@ def config(ctx, mask_cirrus, mask_shadows, mask_method, prob, dark, shadow_dist,
             ctx.obj.cloud_kwargs[key] = val
 
 
-cli.add_command(config)
-
-
 # search command
-@click.command(cls=ChainedCommand)
+@cli.command(cls=ChainedCommand)
 @click.option(
-    '-c', '--collection', type=click.STRING, required=True, callback=_collection_cb,
-    help=f'Earth Engine image collection to search. geedim or EE collection names can be used.'
+    '-c',
+    '--collection',
+    type=click.STRING,
+    required=True,
+    callback=_collection_cb,
+    help=f'Earth Engine image collection to search. geedim or EE collection names can be used.',
 )
 @click.option('-s', '--start-date', type=click.DateTime(), required=False, default=None, help='Start date (UTC).')
 @click.option(
-    '-e', '--end-date', type=click.DateTime(), required=False, default=None,
-    show_default='one day after ``--start-date``', help='End date (UTC).'
+    '-e',
+    '--end-date',
+    type=click.DateTime(),
+    required=False,
+    default=None,
+    show_default='one day after ``--start-date``',
+    help='End date (UTC).',
 )
 @bbox_option
 @click.option(
-    '-r', '--region', type=click.Path(exists=True, dir_okay=False, allow_dash=True), callback=_region_cb,
-    help='Region defined by geojson polygon or raster file. Use "-" to read geojson from stdin.'
-)
-@click.option(
-    '-fp', '--fill-portion', '--fill', metavar='VALUE', type=click.FloatRange(min=0, max=100), default=None,
-    is_flag=False, flag_value=0,
+    '-r',
+    '--region',
+    type=click.Path(exists=True, dir_okay=False, allow_dash=True),
+    callback=_region_cb,
+    help='Region defined by geojson polygon or raster file. Use "-" to read geojson from stdin.',
+)
+@click.option(
+    '-fp',
+    '--fill-portion',
+    '--fill',
+    metavar='VALUE',
+    type=click.FloatRange(min=0, max=100),
+    default=None,
+    is_flag=False,
+    flag_value=0,
     show_default='don\'t calculate, or filter on, fill portion',
     help='Lower limit on the portion of the region that contains filled/valid image pixels (%).  Uses zero if VALUE is '
-         'not specified.'
-)
-@click.option(
-    '-cp', '--cloudless-portion', '--cloudless', metavar='VALUE', type=click.FloatRange(min=0, max=100), default=None,
-    is_flag=False, flag_value=0, show_default='don\'t calculate, or filter on, cloudless portion',
+    'not specified.',
+)
+@click.option(
+    '-cp',
+    '--cloudless-portion',
+    '--cloudless',
+    metavar='VALUE',
+    type=click.FloatRange(min=0, max=100),
+    default=None,
+    is_flag=False,
+    flag_value=0,
+    show_default='don\'t calculate, or filter on, cloudless portion',
     help='Lower limit on the portion of filled pixels that are cloud/shadow free (%).  Uses zero if VALUE is not '
-         'specified.  If cloud/shadow masking is not supported for the specified collection, ``--cloudless-portion`` '
-         'has no effect.'
+    'specified.  If cloud/shadow masking is not supported for the specified collection, ``--cloudless-portion`` '
+    'has no effect.',
 )
 @click.option(
-    '-cf', '--custom-filter', type=click.STRING, default=None,
-    help='Custom image property filter e.g. "property > value".  Quote delimiters are required.'
+    '-cf',
+    '--custom-filter',
+    type=click.STRING,
+    default=None,
+    help='Custom image property filter e.g. "property > value".  Quote delimiters are required.',
 )
 @click.option(
-    '-ap', '--add-property', 'add_props', type=click.STRING, default=None, multiple=True,
-    help='Additional image property name(s) to include in search results.'
+    '-ap',
+    '--add-property',
+    'add_props',
+    type=click.STRING,
+    default=None,
+    multiple=True,
+    help='Additional image property name(s) to include in search results.',
 )
 @click.option(
-    '-op', '--output', type=click.Path(exists=False, dir_okay=False, writable=True), default=None,
-    help='JSON file to write search results to.'
+    '-op',
+    '--output',
+    type=click.Path(exists=False, dir_okay=False, writable=True),
+    default=None,
+    help='JSON file to write search results to.',
 )
 @click.pass_obj
 def search(
-    obj, collection, start_date, end_date, bbox, region, fill_portion, cloudless_portion, custom_filter, output,
-    add_props
+    obj,
+    collection,
+    start_date,
+    end_date,
+    bbox,
+    region,
+    fill_portion,
+    cloudless_portion,
+    custom_filter,
+    output,
+    add_props,
 ):
     # @formatter:off
     """
@@ -488,8 +609,13 @@ def search(
     label = f'Searching for {collection} images: '
     with Spinner(label=label, leave=' '):
         gd_collection = gd_collection.search(
-            start_date, end_date, obj.region, fill_portion=fill_portion, cloudless_portion=cloudless_portion,
-            custom_filter=custom_filter, **obj.cloud_kwargs
+            start_date,
+            end_date,
+            obj.region,
+            fill_portion=fill_portion,
+            cloudless_portion=cloudless_portion,
+            custom_filter=custom_filter,
+            **obj.cloud_kwargs,
         )
         # retrieve search result properties from EE
         num_images = len(gd_collection.properties)
@@ -509,11 +635,8 @@ def search(
             json.dump(gd_collection.properties, f)
 
 
-cli.add_command(search)
-
-
 # download command
-@click.command(cls=ChainedCommand)
+@cli.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
@@ -524,23 +647,41 @@ def search(
 @like_option
 @dtype_option
 @click.option(
-    '-bn', '--band-name', 'bands', type=click.STRING, multiple=True, default=None, show_default='all bands',
-    help='Band name(s) to download.'
+    '-bn',
+    '--band-name',
+    'bands',
+    type=click.STRING,
+    multiple=True,
+    default=None,
+    show_default='all bands',
+    help='Band name(s) to download.',
 )
 @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.'
+    '-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).'
+    '-mts',
+    '--max-tile-size',
+    type=click.FLOAT,
+    default=BaseImage._ee_max_tile_size,
+    show_default=True,
+    help='Maximum download tile size (MB).',
 )
 @click.option(
-    '-mtd', '--max-tile-dim', type=click.INT, default=BaseImage._ee_max_tile_dim, show_default=True,
-    help='Maximum download tile dimension (pixels).'
+    '-mtd',
+    '--max-tile-dim',
+    type=click.INT,
+    default=BaseImage._ee_max_tile_dim,
+    show_default=True,
+    help='Maximum download tile dimension (pixels).',
 )
 @click.option('-o', '--overwrite', is_flag=True, default=False, help='Overwrite the destination file if it exists.')
 @click.pass_obj
@@ -604,26 +745,36 @@ def download(obj, image_id, bbox, region, like, download_dir, mask, max_tile_siz
     for im in image_list:
         filename = pathlib.Path(download_dir).joinpath(im.name + '.tif')
         im.download(
-            filename, region=obj.region, max_tile_size=max_tile_size, max_tile_dim=max_tile_dim, overwrite=overwrite,
-            **kwargs
+            filename,
+            region=obj.region,
+            max_tile_size=max_tile_size,
+            max_tile_dim=max_tile_dim,
+            overwrite=overwrite,
+            **kwargs,
         )
 
 
-cli.add_command(download)
-
-
 # export command
-@click.command(cls=ChainedCommand)
+@cli.command(cls=ChainedCommand)
 @click.option('-i', '--id', 'image_id', type=click.STRING, multiple=True, help='Earth Engine image ID(s) to export.')
 @click.option(
-    '-t', '--type', type=click.Choice([t.value for t in ExportType], case_sensitive=True),
-    default=BaseImage._default_export_type.value, show_default=True, callback=_export_type_cb,
-    help='Export type.'
-)
-@click.option(
-    '-f', '-df', '--folder', '--drive-folder', type=click.STRING, default=None,
+    '-t',
+    '--type',
+    type=click.Choice([t.value for t in ExportType], case_sensitive=True),
+    default=BaseImage._default_export_type.value,
+    show_default=True,
+    callback=_export_type_cb,
+    help='Export type.',
+)
+@click.option(
+    '-f',
+    '-df',
+    '--folder',
+    '--drive-folder',
+    type=click.STRING,
+    default=None,
     help='Google Drive folder, Earth Engine asset project, or Google Cloud Storage bucket to export image(s) to.  '
-         'Interpretation based on :option:`--type`.'
+    'Interpretation based on :option:`--type`.',
 )
 @crs_option
 @bbox_option
@@ -634,8 +785,14 @@ def download(obj, image_id, bbox, region, like, download_dir, mask, max_tile_siz
 @like_option
 @dtype_option
 @click.option(
-    '-bn', '--band-name', 'bands', type=click.STRING, multiple=True, default=None, show_default='all bands',
-    help='Band name(s) to download.'
+    '-bn',
+    '--band-name',
+    'bands',
+    type=click.STRING,
+    multiple=True,
+    default=None,
+    show_default='all bands',
+    help='Band name(s) to download.',
 )
 @mask_option
 @resampling_option
@@ -715,44 +872,63 @@ def export(obj, image_id, type, folder, bbox, region, like, mask, wait, **kwargs
                 obj.image_list += [asset_id(im.name, folder)]
 
 
-cli.add_command(export)
-
-
 # composite command
-@click.command(cls=ChainedCommand)
+@cli.command(cls=ChainedCommand)
 @click.option(
     '-i', '--id', 'image_id', type=click.STRING, multiple=True, help='Earth Engine image ID(s) to include in composite.'
 )
 @click.option(
-    '-cm', '--method', 'method', type=click.Choice([cm.value for cm in CompositeMethod], case_sensitive=False),
-    default=None, callback=_comp_method_cb,
+    '-cm',
+    '--method',
+    'method',
+    type=click.Choice([cm.value for cm in CompositeMethod], case_sensitive=False),
+    default=None,
+    callback=_comp_method_cb,
     show_default='`q-mosaic` for cloud/shadow mask supported collections, `mosaic` otherwise.',
-    help='Compositing method to use.'
+    help='Compositing method to use.',
 )
 @click.option(
-    '-m/-nm', '--mask/--no-mask', default=True, show_default=True,
+    '-m/-nm',
+    '--mask/--no-mask',
+    default=True,
+    show_default=True,
     help='Whether to apply cloud/shadow (or fill) masks to input images before compositing.  Fill masks are used for '
-    'images without support for cloud/shadow masking.'
+    'images without support for cloud/shadow masking.',
 )
 @click.option(
-    '-rs', '--resampling', type=click.Choice([rm.value for rm in ResamplingMethod], case_sensitive=True),
-    default=BaseImage._default_resampling.value, callback=_resampling_method_cb, show_default=True,
-    help='Resample images with this method before compositing.'
+    '-rs',
+    '--resampling',
+    type=click.Choice([rm.value for rm in ResamplingMethod], case_sensitive=True),
+    default=BaseImage._default_resampling.value,
+    callback=_resampling_method_cb,
+    show_default=True,
+    help='Resample images with this method before compositing.',
 )
 @click.option(
-    '-b', '--bbox', type=click.FLOAT, nargs=4, default=None, callback=_bbox_cb,
+    '-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 ``--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,
+    '-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 ``--method``.'
+    'or raster file, region.  Valid for `mosaic` and `q-mosaic` compositing ``--method``.',
 )
 @click.option(
-    '-d', '--date', type=click.DateTime(),
+    '-d',
+    '--date',
+    type=click.DateTime(),
     help='Give preference to images closest to this date (UTC).  Valid for `mosaic` and `q-mosaic` compositing '
-    '``--method``.'
+    '``--method``.',
 )
 @click.pass_obj
 def composite(obj, image_id, mask, method, resampling, bbox, region, date):
@@ -834,6 +1010,5 @@ def composite(obj, image_id, mask, method, resampling, bbox, region, date):
     ]
 
 
-cli.add_command(composite)
-
-##
+if __name__ == '__main__':
+    cli()
diff --git a/geedim/download.py b/geedim/download.py
index 5db248f..c145e1b 100644
--- a/geedim/download.py
+++ b/geedim/download.py
@@ -14,7 +14,6 @@
    limitations under the License.
 """
 
-##
 import logging
 import os
 import pathlib
@@ -29,9 +28,9 @@
 import ee
 import numpy as np
 import rasterio as rio
+from rasterio import windows, features, warp
 from rasterio.crs import CRS
 from rasterio.enums import Resampling as RioResampling
-from rasterio import windows, features, warp
 from tqdm import TqdmWarning
 from tqdm.auto import tqdm
 from tqdm.contrib.logging import logging_redirect_tqdm
@@ -101,26 +100,26 @@ def from_id(cls, image_id: str) -> 'BaseImage':
 
     @property
     def _ee_info(self) -> Dict:
-        """ Earth Engine image metadata. """
+        """Earth Engine image metadata."""
         if self.__ee_info is None:
             self.__ee_info = self._ee_image.getInfo()
         return self.__ee_info
 
     @property
     def _min_projection(self) -> Dict:
-        """ Projection information corresponding to the minimum scale band. """
+        """Projection information corresponding to the minimum scale band."""
         if not self.__min_projection:
             self.__min_projection = self._get_projection(self._ee_info, min_scale=True)
         return self.__min_projection
 
     @property
     def _stac(self) -> Optional[StacItem]:
-        """ Image STAC info.  None if there is no Earth Engine STAC entry for the image / image's collection. """
+        """Image STAC info.  None if there is no Earth Engine STAC entry for the image / image's collection."""
         return StacCatalog().get_item(self.id)
 
     @property
     def ee_image(self) -> ee.Image:
-        """ Encapsulated Earth Engine image. """
+        """Encapsulated Earth Engine image."""
         return self._ee_image
 
     @ee_image.setter
@@ -132,7 +131,7 @@ def ee_image(self, value: ee.Image):
 
     @property
     def id(self) -> Optional[str]:
-        """ Earth Engine image ID.  None if the image ``system:id`` property is not set. """
+        """Earth Engine image ID.  None if the image ``system:id`` property is not set."""
         if self._id:  # avoid a call to getInfo() if _id is set
             return self._id
         else:
@@ -140,7 +139,7 @@ def id(self) -> Optional[str]:
 
     @property
     def date(self) -> Optional[datetime]:
-        """ Image capture date & time.  None if the image ``system:time_start`` property is not set. """
+        """Image capture date & time.  None if the image ``system:time_start`` property is not set."""
         if 'system:time_start' in self.properties:
             return datetime.utcfromtimestamp(self.properties['system:time_start'] / 1000)
         else:
@@ -148,7 +147,7 @@ def date(self) -> Optional[datetime]:
 
     @property
     def name(self) -> Optional[str]:
-        """ Image name (the :attr:`id` with slashes replaced by dashes). """
+        """Image name (the :attr:`id` with slashes replaced by dashes)."""
         return self.id.replace('/', '-') if self.id else None
 
     @property
@@ -160,12 +159,12 @@ def crs(self) -> Optional[str]:
 
     @property
     def scale(self) -> Optional[float]:
-        """ Minimum scale (m) of image bands. None if the image has no fixed projection. """
+        """Minimum scale (m) of image bands. None if the image has no fixed projection."""
         return self._min_projection['scale']
 
     @property
     def footprint(self) -> Optional[Dict]:
-        """ Geojson polygon of the image extent.  None if the image is a composite. """
+        """Geojson polygon of the image extent.  None if the image is a composite."""
         if ('properties' not in self._ee_info) or ('system:footprint' not in self._ee_info['properties']):
             return None
         footprint = self._ee_info['properties']['system:footprint']
@@ -181,29 +180,29 @@ def footprint(self) -> Optional[Dict]:
 
     @property
     def shape(self) -> Optional[Tuple[int, int]]:
-        """ (row, column) pixel dimensions of the minimum scale band. None if the image has no fixed projection. """
+        """(row, column) pixel dimensions of the minimum scale band. None if the image has no fixed projection."""
         return self._min_projection['shape']
 
     @property
     def count(self) -> int:
-        """ Number of image bands. """
+        """Number of image bands."""
         return len(self._ee_info['bands']) if 'bands' in self._ee_info else None
 
     @property
     def transform(self) -> Optional[rio.Affine]:
-        """ Geo-transform of the minimum scale band. None if the image has no fixed projection. """
+        """Geo-transform of the minimum scale band. None if the image has no fixed projection."""
         return self._min_projection['transform']
 
     @property
     def dtype(self) -> str:
-        """ Minimum size data type able to represent the image values. """
+        """Minimum size data type able to represent the image values."""
         if not self._min_dtype:
             self._min_dtype = self._get_min_dtype(self._ee_info)
         return self._min_dtype
 
     @property
     def size(self) -> Optional[int]:
-        """ Image size (bytes).  None if the image has no fixed projection. """
+        """Image size (bytes).  None if the image has no fixed projection."""
         if not self.shape:
             return None
         dtype_size = np.dtype(self.dtype).itemsize
@@ -211,37 +210,41 @@ def size(self) -> Optional[int]:
 
     @property
     def has_fixed_projection(self) -> bool:
-        """ True if the image has a fixed projection, otherwise False. """
+        """True if the image has a fixed projection, otherwise False."""
         return self.scale is not None
 
     @property
     def properties(self) -> Dict:
-        """ Earth Engine image properties. """
+        """Earth Engine image properties."""
         return self._ee_info['properties'] if 'properties' in self._ee_info else {}
 
     @property
     def band_properties(self) -> List[Dict]:
-        """ Merged STAC and Earth Engine band properties. """
+        """Merged STAC and Earth Engine band properties."""
         return self._get_band_properties()
 
     @property
     def refl_bands(self) -> Optional[List[str]]:
-        """ List of spectral / reflectance bands, if any. """
+        """List of spectral / reflectance bands, if any."""
         if not self._stac:
             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. """
+        """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
+            crs=self.crs,
+            transform=self.transform,
+            width=self.shape[1],
+            height=self.shape[0],
+            dtype=self.dtype,
+            count=self.count,
         )
 
     @property
     def bounded(self) -> bool:
-        """ True if the image is bounded, otherwise False. """
+        """True if the image is bounded, otherwise False."""
         # TODO: an unbounded region could also have these bounds
         unbounded_bounds = (-180, -90, 180, 90)
         return (self.footprint is not None) and (features.bounds(self.footprint) != unbounded_bounds)
@@ -272,12 +275,15 @@ def _get_projection(ee_info: Dict, min_scale=True) -> Dict:
 
     @staticmethod
     def _get_min_dtype(ee_info: Dict) -> str:
-        """ Return the minimum size rasterio data type corresponding to a given Earth Engine image info dictionary. """
+        """Return the minimum size rasterio data type corresponding to a given Earth Engine image info dictionary."""
+
         def get_min_int_dtype(band_info: List[Dict]) -> Optional[str]:
-            """ Return the minimum integer dtype for the given band information. """
+            """Return the minimum integer dtype for the given band information."""
             # create a list of integer min/max values
             int_minmax = [
-                minmax for band_dict in band_info if band_dict['data_type']['precision'] == 'int'
+                minmax
+                for band_dict in band_info
+                if band_dict['data_type']['precision'] == 'int'
                 for minmax in (int(band_dict['data_type']['min']), int(band_dict['data_type']['max']))
             ]  # yapf: disable
 
@@ -319,10 +325,16 @@ def _str_format_size(byte_size: float, units=['bytes', 'KB', 'MB', 'GB', 'TB', '
 
     @staticmethod
     def _convert_dtype(ee_image: ee.Image, dtype: str) -> ee.Image:
-        """ Convert the data type of an Earth Engine image to a specified type. """
+        """Convert the data type of an Earth Engine image to a specified type."""
         conv_dict = dict(
-            float32=ee.Image.toFloat, float64=ee.Image.toDouble, uint8=ee.Image.toUint8, int8=ee.Image.toInt8,
-            uint16=ee.Image.toUint16, int16=ee.Image.toInt16, uint32=ee.Image.toUint32, int32=ee.Image.toInt32,
+            float32=ee.Image.toFloat,
+            float64=ee.Image.toDouble,
+            uint8=ee.Image.toUint8,
+            int8=ee.Image.toInt8,
+            uint16=ee.Image.toUint16,
+            int16=ee.Image.toInt16,
+            uint32=ee.Image.toUint32,
+            int32=ee.Image.toInt32,
         )
         if dtype not in conv_dict:
             raise TypeError(f'Unsupported dtype: {dtype}')
@@ -330,7 +342,7 @@ def _convert_dtype(ee_image: ee.Image, dtype: str) -> ee.Image:
         return conv_dict[dtype](ee_image)
 
     def _get_band_properties(self) -> List[Dict]:
-        """ Merge Earth Engine and STAC band properties for this image. """
+        """Merge Earth Engine and STAC band properties for this image."""
         band_ids = [bd['id'] for bd in self._ee_info['bands']]
         if self._stac:
             stac_bands_props = self._stac.band_props
@@ -365,8 +377,8 @@ def _scale_offset(ee_image: ee.Image, band_properties: List[Dict]) -> ee.Image:
             return ee_image
 
         # make band scale and offset dicts
-        scale_dict = {bp['name']: bp['scale'] if 'scale' in bp else 1. for bp in band_properties}
-        offset_dict = {bp['name']: bp['offset'] if 'offset' in bp else 0. for bp in band_properties}
+        scale_dict = {bp['name']: bp['scale'] if 'scale' in bp else 1.0 for bp in band_properties}
+        offset_dict = {bp['name']: bp['offset'] if 'offset' in bp else 0.0 for bp in band_properties}
 
         if all([s == 1 for s in scale_dict.values()]) and all([o == 0 for o in offset_dict.values()]):
             # all scales==1 and all offsets==0
@@ -389,9 +401,16 @@ 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, 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, bands: List[str] = None,
+        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,
+        bands: List[str] = None,
     ) -> 'BaseImage':
         """
         Prepare the encapsulated image for export/download.  Will reproject, resample, clip and convert the image
@@ -452,8 +471,7 @@ def _prepare_for_export(
             # 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.'
+                f'This image is unbounded, you need to specify a region; or a crs, crs_transform and shape.'
             )
 
         if not scale and not shape and self.crs == 'EPSG:4326':
@@ -527,8 +545,14 @@ def _prepare_for_export(
         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, bands=bands,
-            fileFormat='GeoTIFF', filePerBand=False
+            crs=crs,
+            crs_transform=crs_transform,
+            dimensions=dimensions,
+            region=region,
+            scale=scale,
+            bands=bands,
+            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}
@@ -559,9 +583,18 @@ def _prepare_for_download(self, set_nodata: bool = True, **kwargs) -> Tuple['Bas
         )  # yapf: disable
         nodata = nodata_dict[exp_image.dtype] if set_nodata else None
         profile = dict(
-            driver='GTiff', dtype=exp_image.dtype, nodata=nodata, width=exp_image.shape[1], height=exp_image.shape[0],
-            count=exp_image.count, crs=CRS.from_string(utils.rio_crs(exp_image.crs)), transform=exp_image.transform,
-            compress='deflate', interleave='band', tiled=True, photometric='MINISBLACK',
+            driver='GTiff',
+            dtype=exp_image.dtype,
+            nodata=nodata,
+            width=exp_image.shape[1],
+            height=exp_image.shape[0],
+            count=exp_image.count,
+            crs=CRS.from_string(utils.rio_crs(exp_image.crs)),
+            transform=exp_image.transform,
+            compress='deflate',
+            interleave='band',
+            tiled=True,
+            photometric='MINISBLACK',
         )
         # add BIGTIFF support if the uncompressed image is bigger than 4GB
         if exp_image.size >= 4e9:
@@ -587,7 +620,7 @@ def _get_tile_shape(
                 f'`max_tile_dim` must be less than or equal to the Earth Engine download limit of '
                 f'{BaseImage._ee_max_tile_size} pixels.'
             )
-        max_tile_dim = max_tile_dim or BaseImage._ee_max_tile_dim   # set max_tile_dim to EE default if None
+        max_tile_dim = max_tile_dim or BaseImage._ee_max_tile_dim  # set max_tile_dim to EE default if None
 
         # find the total number of tiles the image must be divided into to satisfy max_tile_size
         image_shape = np.array(self.shape, dtype='int64')
@@ -615,7 +648,7 @@ def _get_tile_shape(
         return tile_shape, num_tiles
 
     def _build_overviews(self, filename: Union[str, pathlib.Path], max_num_levels: int = 8, min_ovw_pixels: int = 256):
-        """ Build internal overviews, downsampled by successive powers of 2, for a given filename. """
+        """Build internal overviews, downsampled by successive powers of 2, for a given filename."""
 
         # TODO: revisit multi-threaded overviews on gdal update
         # build overviews in a single threaded environment (currently gdal reports errors when building overviews
@@ -635,12 +668,12 @@ def _build_overviews(self, filename: Union[str, pathlib.Path], max_num_levels: i
                 ds.build_overviews(ovw_levels, RioResampling.average)
 
     def _write_metadata(self, dataset: rio.io.DatasetWriter):
-        """ Write Earth Engine and STAC metadata to an open rasterio dataset. """
+        """Write Earth Engine and STAC metadata to an open rasterio dataset."""
         if dataset.closed:
             raise IOError('Image dataset is closed')
 
         # replace 'system:*' property keys with 'system-*', and remove footprint if its there
-        properties = {k.replace(':', '-'):v for k, v in self.properties.items()}
+        properties = {k.replace(':', '-'): v for k, v in self.properties.items()}
         if 'system-footprint' in properties:
             properties.pop('system-footprint')
         dataset.update_tags(**properties)
@@ -649,7 +682,7 @@ def _write_metadata(self, dataset: rio.io.DatasetWriter):
             dataset.update_tags(LICENSE=self._stac.license)
 
         def clean_text(text, width=80) -> str:
-            """ Return a shortened tidied string. """
+            """Return a shortened tidied string."""
             if not isinstance(text, str):
                 return text
             text = text.split('.')[0] if len(text) > width else text
@@ -807,8 +840,12 @@ def export(
         type = ExportType(type)
         if type == ExportType.drive:
             task = ee.batch.Export.image.toDrive(
-                image=exp_image.ee_image, description=filename[:100], folder=folder, fileNamePrefix=filename,
-                maxPixels=1e9, formatOptions=dict(cloudOptimized=True),
+                image=exp_image.ee_image,
+                description=filename[:100],
+                folder=folder,
+                fileNamePrefix=filename,
+                maxPixels=1e9,
+                formatOptions=dict(cloudOptimized=True),
             )
         elif type == ExportType.asset:
             # if folder is specified create an EE asset ID from it and filename,
@@ -817,11 +854,14 @@ def export(
             # fix description for when filename is asset id with forward slashes
             description = filename.replace('/', '-')[:100]
             task = ee.batch.Export.image.toAsset(
-                image=exp_image.ee_image, description=description, assetId=asset_id, maxPixels=1e9,
+                image=exp_image.ee_image, description=description, assetId=asset_id, maxPixels=1e9
             )
         else:
             task = ee.batch.Export.image.toCloudStorage(
-                image=exp_image.ee_image, description=filename[:100], bucket=folder, maxPixels=1e9,
+                image=exp_image.ee_image,
+                description=filename[:100],
+                bucket=folder,
+                maxPixels=1e9,
                 formatOptions=dict(cloudOptimized=True),
             )
 
@@ -831,8 +871,13 @@ def export(
         return task
 
     def download(
-        self, filename: Union[pathlib.Path, str], overwrite: bool = False, num_threads: Optional[int] = None,
-        max_tile_size: Optional[float] = None, max_tile_dim: Optional[int] = None, **kwargs
+        self,
+        filename: Union[pathlib.Path, str],
+        overwrite: bool = False,
+        num_threads: Optional[int] = None,
+        max_tile_size: Optional[float] = None,
+        max_tile_dim: Optional[int] = None,
+        **kwargs,
     ):
         """
         Download the encapsulated image to a GeoTiff file.
@@ -930,7 +975,7 @@ def download(
             desc=desc, total=raw_download_size, bar_format=bar_format, dynamic_ncols=True, unit_scale=True, unit='B'
         )
 
-        session = utils.retry_session(5)
+        session = utils.retry_session()
         warnings.filterwarnings('ignore', category=TqdmWarning)
         # redirect logging through tqdm
         redir_tqdm = logging_redirect_tqdm([logging.getLogger(__package__)], tqdm_class=type(bar))
@@ -938,7 +983,7 @@ def download(
         with redir_tqdm, env, rio.open(filename, 'w', **profile) as out_ds, bar:
 
             def download_tile(tile):
-                """Download a tile and write into the destination GeoTIFF. """
+                """Download a tile and write into the destination GeoTIFF."""
                 tile_array = tile.download(session=session, bar=bar)
                 with out_lock:
                     out_ds.write(tile_array, window=tile.window)
@@ -955,7 +1000,7 @@ def download_tile(tile):
                     executor.shutdown(wait=False, cancel_futures=True)
                     raise ex
 
-            bar.update(bar.total - bar.n)   # ensure the bar reaches 100%
+            bar.update(bar.total - bar.n)  # ensure the bar reaches 100%
             # populate GeoTIFF metadata
             exp_image._write_metadata(out_ds)
 
diff --git a/geedim/errors.py b/geedim/errors.py
index 21206b6..36c52db 100644
--- a/geedim/errors.py
+++ b/geedim/errors.py
@@ -16,12 +16,16 @@
 
 
 class GeedimError(Exception):
-    """ Base exception class. """
+    """Base exception class."""
 
 
 class UnfilteredError(GeedimError):
-    """ Raised when attempting to retrieve the properties of an unfiltered image collection. """
+    """Raised when attempting to retrieve the properties of an unfiltered image collection."""
 
 
 class InputImageError(GeedimError):
-    """ Raised when there is a problem with the images making up a collection. """
+    """Raised when there is a problem with the images making up a collection."""
+
+
+class TileError(GeedimError):
+    """Raised when there is a problem downloading an image tile."""
diff --git a/geedim/tile.py b/geedim/tile.py
index 15b8d7a..d683ad1 100644
--- a/geedim/tile.py
+++ b/geedim/tile.py
@@ -14,17 +14,24 @@
    limitations under the License.
 """
 
+import logging
+import threading
+import time
 import zipfile
 from io import BytesIO
-import threading
 
 import numpy as np
-import requests
 import rasterio as rio
+import requests
 from rasterio import Affine, MemoryFile
 from rasterio.windows import Window
+from requests.exceptions import RequestException, JSONDecodeError
 from tqdm.auto import tqdm
+
 from geedim import utils
+from geedim.errors import TileError
+
+logger = logging.getLogger(__name__)
 
 
 class Tile:
@@ -51,22 +58,89 @@ def __init__(self, exp_image, window: Window):
 
     @property
     def window(self) -> Window:
-        """ rasterio tile window into the source image. """
+        """rasterio tile window into the source image."""
         return self._window
 
-    def _get_download_url_response(self, session=None):
-        """ Get tile download url and response. """
-        session = session if session else requests
-        with self._ee_lock:
-            url = self._exp_image.ee_image.getDownloadURL(
-                dict(
-                    crs=self._exp_image.crs, crs_transform=tuple(self._transform)[:6], dimensions=self._shape[::-1],
-                    filePerBand=False, fileFormat='GeoTIFF'
-                )
-            )
-        return session.get(url, stream=True, timeout=(30, 300)), url
+    @staticmethod
+    def _raise_for_status(response: requests.Response):
+        """Raise a TileError if the tile cannot be downloaded."""
+        download_size = int(response.headers.get('content-length', 0))
+        if download_size == 0 or not response.status_code == 200:
+            msg = f'Error downloading tile: {response.status_code} - {response.reason}. URL: {response.url}.'
+            try:
+                resp_dict = response.json()
+                if 'error' in resp_dict and 'message' in resp_dict['error']:
+                    # raise an exception with the response error message
+                    msg = resp_dict['error']['message']
+                    msg = f'Error downloading tile: {msg} URL: {response.url}.'
+                    if 'user memory limit exceeded' in msg.lower():
+                        msg += (
+                            '\nThe `max_tile_size` or `max_tile_dim` parameters can be decreased to work around this '
+                            'error.  Alternatively you can export to Earth Engine asset, and then download the asset '
+                            'image.'
+                        )
+            except JSONDecodeError:
+                pass
+
+            raise TileError(msg)
+
+    def _download_to_array(self, url: str, session: requests.Session = None, bar: tqdm = None) -> np.ndarray:
+        """Download the image tile into a numpy array."""
+        # get image download response
+        session = session or requests
+        response = session.get(url, stream=True, timeout=(30, 300))
+
+        # raise a TileError if the tile cannot be downloaded
+        self._raise_for_status(response)
+
+        # find raw download size
+        download_size = int(response.headers.get('content-length', 0))
+        dtype_size = np.dtype(self._exp_image.dtype).itemsize
+        raw_download_size = self._shape[0] * self._shape[1] * self._exp_image.count * dtype_size
+
+        # download and unzip the tile
+        downloaded_size = 0
+        try:
+            # download zip into buffer
+            zip_buffer = BytesIO()
+            for data in response.iter_content(chunk_size=10240):
+                if bar:
+                    # update with raw download progress (0-1)
+                    bar.update(raw_download_size * (len(data) / download_size))
+                zip_buffer.write(data)
+                downloaded_size += len(data)
+            zip_buffer.flush()
+
+            # extract GeoTIFF bytes from zipped buffer
+            with zipfile.ZipFile(zip_buffer, 'r') as zip_file:
+                gtiff_bytes = zip_file.read(zip_file.filelist[0])
+
+        except (RequestException, zipfile.BadZipfile):
+            if bar:
+                # reverse progress bar
+                bar.update(-raw_download_size * (downloaded_size / download_size))
+                pass
+            raise
+
+        # read the tile array from the GeoTIFF bytes, via a rasterio memory file
+        env = rio.Env(GDAL_NUM_THREADS='ALL_CPUs', GTIFF_FORCE_RGBA=False)
+        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(self, session=None, response=None, bar: tqdm = None):
+    def download(
+        self,
+        session: requests.Session = None,
+        bar: tqdm = None,
+        max_retries: int = 5,
+        backoff_factor: float = 2.0,
+    ) -> np.ndarray:
         """
         Download the image tile into a numpy array.
 
@@ -74,74 +148,50 @@ def download(self, session=None, response=None, bar: tqdm = None):
         ----------
         session: requests.Session, optional
             requests session to use for downloading
-        response: requests.Response, optional
-            Response to a get request on the tile download url.
         bar: tqdm, optional
             tqdm progress bar instance to update with incremental (0-1) download progress.
+        max_retries: int, optional
+            Number of times to retry downloading the tile.  This is independent of the ``session``, which may have its
+            own retry configuration.
+        backoff_factor: float, optional
+            Backoff factor to apply between tile download retries.  The delay between retries is: {backoff_factor} *
+            (2 ** ({number of previous retries})) seconds.  This is independent of the ``session``, which may have its
+            own retry configuration.
 
         Returns
         -------
         array: numpy.ndarray
             3D numpy array of the tile pixel data with bands down the first dimension.
         """
+        session = session or requests
 
-        # get image download url and response
-        if response is None:
-            response, url = self._get_download_url_response(session=session)
-
-        # find raw and actual download sizes
-        dtype_size = np.dtype(self._exp_image.dtype).itemsize
-        raw_download_size = self._shape[0] * self._shape[1] * self._exp_image.count * dtype_size
-        download_size = int(response.headers.get('content-length', 0))
-
-        if download_size == 0 or not response.status_code == 200:
-            resp_dict = response.json()
-            if 'error' in resp_dict and 'message' in resp_dict['error']:
-                msg = resp_dict['error']['message']
-                ex_msg = f'Error downloading tile: {msg}'
-                if 'user memory limit exceeded' in msg.lower():
-                    ex_msg += (
-                        '\nThe `max_tile_size` or `max_tile_dim` parameters can be decreased to work around this '
-                        'error.  Alternatively you can export to Earth Engine asset, and then download the asset image.'
-                    )
-            else:
-                ex_msg = str(response.json())
-            raise IOError(ex_msg)
-
-        # download zip into buffer
-        zip_buffer = BytesIO()
-        downloaded_size = 0
-        with response:
-            for data in response.iter_content(chunk_size=10240):
-                zip_buffer.write(data)
-                if bar is not None:
-                    # update with raw download progress (0-1)
-                    bar.update(raw_download_size * (len(data) / download_size))
-                downloaded_size += len(data)
-            zip_buffer.flush()
-
-        # check downloaded size matches the header content-length
-        if downloaded_size < download_size:
-            raise IOError(
-                f'Incomplete read: {downloaded_size} bytes read, {download_size-downloaded_size} '
-                f'more expected. URL: {response.url}.'
+        # get download URL
+        with self._ee_lock:
+            url = self._exp_image.ee_image.getDownloadURL(
+                dict(
+                    crs=self._exp_image.crs,
+                    crs_transform=tuple(self._transform)[:6],
+                    dimensions=self._shape[::-1],
+                    filePerBand=False,
+                    fileFormat='GeoTIFF',
+                )
             )
 
-        # extract geotiff from zipped buffer into another buffer
-        try:
-            zip_file = zipfile.ZipFile(zip_buffer)
-            ext_buffer = BytesIO(zip_file.read(zip_file.filelist[0]))
-        except zipfile.BadZipfile as ex:
-            raise IOError(f'Bad tile zip from URL: {response.url}.') from ex
-
-        # read the geotiff with a rasterio memory file
-        env = rio.Env(GDAL_NUM_THREADS='ALL_CPUs', GTIFF_FORCE_RGBA=False)
-        with utils.suppress_rio_logs(), env, MemoryFile(ext_buffer) 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
+        # download and read the tile, with retries
+        array = None
+        for retry in range(max_retries + 1):
+            try:
+                array = self._download_to_array(url, session=session, bar=bar)
+                break
+
+            except (RequestException, zipfile.BadZipfile) as ex:
+                if retry < max_retries:
+                    # retry tile download
+                    time.sleep(backoff_factor * (2**retry))
+                    logger.warning(
+                        f'Tile downloaded failed, retry {retry + 1} of {max_retries}.  URL: {url}. {str(ex)}.'
+                    )
+                else:
+                    raise TileError(f'Tile download failed, reached the maximum retries.  URL: {url}.') from ex
 
         return array
diff --git a/geedim/utils.py b/geedim/utils.py
index d8b2e07..c5f00c9 100644
--- a/geedim/utils.py
+++ b/geedim/utils.py
@@ -13,28 +13,30 @@
    See the License for the specific language governing permissions and
    limitations under the License.
 """
-import logging
-from contextlib import contextmanager
+
 import itertools
 import json
+import logging
 import os
 import pathlib
 import sys
 import time
+from contextlib import contextmanager
 from threading import Thread
 from typing import Tuple, Optional
 
 import ee
+import numpy as np
 import rasterio as rio
 import requests
-from geedim.enums import ResamplingMethod
-from rasterio.warp import transform_geom
 from rasterio.env import GDALVersion
+from rasterio.warp import transform_geom
 from rasterio.windows import Window
 from requests.adapters import HTTPAdapter, Retry
-import numpy as np
 from tqdm.auto import tqdm
 
+from geedim.enums import ResamplingMethod
+
 if '__file__' in globals():
     root_path = pathlib.Path(__file__).absolute().parents[1]
 else:
@@ -42,7 +44,8 @@
 
 _GDAL_AT_LEAST_35 = GDALVersion.runtime().at_least("3.5")
 
-def Initialize(opt_url: Optional[str]='https://earthengine-highvolume.googleapis.com', **kwargs):
+
+def Initialize(opt_url: Optional[str] = 'https://earthengine-highvolume.googleapis.com', **kwargs):
     """
     Initialise Earth Engine.
 
@@ -76,7 +79,7 @@ def Initialize(opt_url: Optional[str]='https://earthengine-highvolume.googleapis
 
 
 def singleton(cls):
-    """ Class decorator to make it a singleton. """
+    """Class decorator to make it a singleton."""
     instances = {}
 
     def getinstance() -> cls:
@@ -110,7 +113,7 @@ def split_id(image_id: str) -> Tuple[str, str]:
 
 @contextmanager
 def suppress_rio_logs(level: int = logging.ERROR):
-    """ A context manager that sets the `rasterio` logging level, then returns it to its original value. """
+    """A context manager that sets the `rasterio` logging level, then returns it to its original value."""
     try:
         # GEE sets GeoTIFF `colorinterp` tags incorrectly. This suppresses `rasterio` warning relating to this:
         # 'Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel'
@@ -149,8 +152,10 @@ def get_bounds(filename: pathlib.Path, expand: float = 5):
             bbox_expand = bbox
 
         coordinates = [
-            [bbox_expand.right, bbox_expand.bottom], [bbox_expand.right, bbox_expand.top],
-            [bbox_expand.left, bbox_expand.top], [bbox_expand.left, bbox_expand.bottom],
+            [bbox_expand.right, bbox_expand.bottom],
+            [bbox_expand.right, bbox_expand.top],
+            [bbox_expand.left, bbox_expand.top],
+            [bbox_expand.left, bbox_expand.bottom],
             [bbox_expand.right, bbox_expand.bottom],
         ]
 
@@ -236,7 +241,7 @@ def __exit__(self, exc_type, exc_val, exc_tb):
         self.join()
 
     def run(self):
-        """ Run the spinner thread. """
+        """Run the spinner thread."""
         cursors_it = itertools.cycle('/-\|')
 
         while self._run:
@@ -254,12 +259,12 @@ def run(self):
         self._file.flush()
 
     def start(self):
-        """ Start the spinner thread. """
+        """Start the spinner thread."""
         self._run = True
         Thread.start(self)
 
     def stop(self):
-        """ Stop the spinner thread. """
+        """Stop the spinner thread."""
         self._run = False
 
 
@@ -292,7 +297,7 @@ def resample(ee_image: ee.Image, method: ResamplingMethod) -> ee.Image:
     has_fixed_proj = proj.crs().compareTo('EPSG:4326').neq(0).Or(proj.nominalScale().toInt64().neq(111319))
 
     def _resample(ee_image: ee.Image) -> ee.Image:
-        """ Resample the given image, allowing for additional 'average' method. """
+        """Resample the given image, allowing for additional 'average' method."""
         if method == ResamplingMethod.average:
             return ee_image.reduceResolution(reducer=ee.Reducer.mean(), maxPixels=1024)
         else:
@@ -302,10 +307,12 @@ def _resample(ee_image: ee.Image) -> ee.Image:
 
 
 def retry_session(
-    retries: int = 5, backoff_factor: float = 2., status_forcelist: Tuple = (429, 500, 502, 503, 504),
-    session: requests.Session = None
+    retries: int = 5,
+    backoff_factor: float = 2.0,
+    status_forcelist: Tuple = (429, 500, 502, 503, 504),
+    session: requests.Session = None,
 ) -> requests.Session:
-    """ requests session configured for retries. """
+    """requests session configured for retries."""
     session = session or requests.Session()
     retry = Retry(
         total=retries, read=retries, connect=retries, backoff_factor=backoff_factor, status_forcelist=status_forcelist
@@ -342,7 +349,7 @@ def expand_window_to_grid(win: Window, expand_pixels: Tuple[int, int] = (0, 0))
 
 
 def rio_crs(crs: str) -> str:
-    """ Convert a GEE CRS string to a rasterio compatible CRS string. """
+    """Convert a GEE CRS string to a rasterio compatible CRS string."""
     if crs == 'SR-ORG:6974':
         # This is a workaround for https://issuetracker.google.com/issues/194561313, that replaces the alleged GEE
         # SR-ORG:6974 with actual (rasterio) SR-ORG:6842. WKT taken from
diff --git a/geedim/version.py b/geedim/version.py
index 6da2111..b280975 100644
--- a/geedim/version.py
+++ b/geedim/version.py
@@ -1,1 +1,1 @@
-__version__ = '1.7.3'
+__version__ = '1.8.0'
diff --git a/pyproject.toml b/pyproject.toml
index b5a3c46..aa48457 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,55 @@
-[build-system]
-requires = [
-    "setuptools>=42",
-    "wheel"
+[project]
+name = 'geedim'
+description = 'Search, composite and download Google Earth Engine imagery.'
+readme = 'README.rst'
+requires-python = '>=3.6'
+dependencies = [
+    'numpy>=1.19',
+    'rasterio>=1.1',
+    'click>=8',
+    'tqdm>=4.6',
+    'earthengine-api>=0.1.2',
+    'requests>=2.2',
+    'tabulate>=0.8'
+]
+authors = [{name = 'Leftfield Geospatial'}]
+keywords = [
+    'earth engine', 'satellite imagery', 'search', 'download', 'composite', 'cloud', 'shadow'
 ]
-build-backend = "setuptools.build_meta"
\ No newline at end of file
+classifiers = [
+    'Programming Language :: Python :: 3',
+    'License :: OSI Approved :: Apache Software License',
+    'Operating System :: OS Independent',
+]
+dynamic = ['version']
+
+[project.scripts]
+geedim = 'geedim.cli:cli'
+
+[project.optional-dependencies]
+tests = ['pytest', 'pytest-xdist']
+docs = ['sphinx', 'sphinx-click', 'sphinx-rtd-theme', 'nbsphinx', 'ipykernel', 'jupyter-sphinx']
+
+[project.urls]
+Homepage = 'https://github.com/leftfield-geospatial/geedim'
+Documentation = "https://geedim.readthedocs.org"
+Source = 'https://github.com/leftfield-geospatial/geedim'
+Changelog = 'https://github.com/leftfield-geospatial/geedim/releases'
+Issues = 'https://github.com/leftfield-geospatial/geedim/issues'
+
+[build-system]
+requires = ['setuptools']
+build-backend = 'setuptools.build_meta'
+
+[tool.setuptools]
+packages = ['geedim']
+
+[tool.setuptools.dynamic]
+version = {attr = 'geedim.version.__version__'}
+
+[tool.setuptools.package-data]
+geedim = ['data/ee_stac_urls.json']
+
+[tool.black]
+line-length = 120
+skip-string-normalization = true
diff --git a/setup.py b/setup.py
deleted file mode 100644
index f0a4353..0000000
--- a/setup.py
+++ /dev/null
@@ -1,74 +0,0 @@
-"""
-   Copyright 2021 Dugal Harris - dugalh@gmail.com
-
-   Licensed under the Apache License, Version 2.0 (the "License");
-   you may not use this file except in compliance with the License.
-   You may obtain a copy of the License at
-
-       http://www.apache.org/licenses/LICENSE-2.0
-
-   Unless required by applicable law or agreed to in writing, software
-   distributed under the License is distributed on an "AS IS" BASIS,
-   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-   See the License for the specific language governing permissions and
-   limitations under the License.
-"""
-
-import sys
-from pathlib import Path
-
-from setuptools import setup, find_packages
-"""
- Build and upload to testpypi:
-     conda install -c conda-forge build twine
-     python -m build
-     python -m twine upload --repository testpypi dist/*
-
- Install from testpypi:
-    python -m pip install --extra-index-url https://test.pypi.org/simple/ geedim
-
- Install local development version:
-    pip install -e .
-"""
-
-this_directory = Path(__file__).parent
-long_description = (this_directory / 'README.rst').read_text()
-sys.path[0:0] = ['geedim']
-from version import __version__
-
-setup(
-    name='geedim',
-    version=__version__,
-    description='Search, composite and download Google Earth Engine imagery.',
-    long_description=long_description,
-    long_description_content_type='text/x-rst',
-    author='Dugal Harris',
-    author_email='dugalh@gmail.com',
-    url='https://github.com/leftfield-geospatial/geedim',
-    license='Apache-2.0',
-    packages=find_packages(include=['geedim']),
-    package_data={'geedim': ['data/ee_stac_urls.json']},
-    install_requires=[
-        'numpy>=1.19',
-        'rasterio>=1.1',
-        'click>=8',
-        'tqdm>=4.6',
-        'earthengine-api>=0.1.2',
-        'requests>=2.2',
-        'tabulate>=0.8',
-    ],
-    python_requires='>=3.6',
-    classifiers=[
-        'Programming Language :: Python :: 3',
-        'License :: OSI Approved :: Apache Software License',
-        'Operating System :: OS Independent',
-    ],
-    keywords=[
-        'earth engine', 'satellite imagery', 'search', 'download', 'composite', 'cloud', 'shadow',
-    ],
-    entry_points={'console_scripts': ['geedim=geedim.cli:cli']},
-    project_urls={
-        'Documentation': 'https://geedim.readthedocs.io',
-        'Source': 'https://github.com/leftfield-geospatial/geedim',
-    },
-)  # yapf: disable


Test patch:
diff --git a/.github/workflows/install-test-conda-forge.yml b/.github/workflows/install-test-conda-forge.yml
index 71022db..9bdbf88 100644
--- a/.github/workflows/install-test-conda-forge.yml
+++ b/.github/workflows/install-test-conda-forge.yml
@@ -29,7 +29,7 @@ jobs:
     - name: Install package
       run: |
         conda info
-        conda install geedim>=1.7.3
+        conda install geedim>=1.8.0
         conda list
 
     - name: Run tests
diff --git a/tests/test_download.py b/tests/test_download.py
index d14277f..6516212 100644
--- a/tests/test_download.py
+++ b/tests/test_download.py
@@ -13,25 +13,27 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
+
+import copy
 import pathlib
 from datetime import datetime
 from typing import Dict, Tuple, List
-import copy
 
 import ee
 import numpy as np
 import pytest
 import rasterio as rio
+from rasterio import Affine, windows
 from rasterio import features, warp
 from rasterio.crs import CRS
+
+from geedim import utils
 from geedim.download import BaseImage
 from geedim.enums import ResamplingMethod, ExportType
-from geedim import utils
-from rasterio import Affine, windows
 
 
 class BaseImageLike:
-    """ Emulate BaseImage for _get_tile_shape() and tiles(). """
+    """Emulate BaseImage for _get_tile_shape() and tiles()."""
 
     def __init__(
         self, shape: Tuple[int, int], count: int = 10, dtype: str = 'uint16', transform: Affine = Affine.identity()
@@ -49,7 +51,7 @@ def __init__(
 
 @pytest.fixture(scope='session')
 def user_base_image() -> BaseImage:
-    """ A BaseImage instance where the encapsulated image has no fixed projection or ID.  """
+    """A BaseImage instance where the encapsulated image has no fixed projection or ID."""
     return BaseImage(ee.Image([1, 2, 3]))
 
 
@@ -64,48 +66,51 @@ def user_fix_base_image() -> BaseImage:
 
 @pytest.fixture(scope='session')
 def s2_sr_base_image(s2_sr_image_id: str) -> BaseImage:
-    """ A BaseImage instance encapsulating a Sentinel-2 image.  Covers `region_*ha`.  """
+    """A BaseImage instance encapsulating a Sentinel-2 image.  Covers `region_*ha`."""
     return BaseImage.from_id(s2_sr_image_id)
 
 
 @pytest.fixture(scope='session')
 def l9_base_image(l9_image_id: str) -> BaseImage:
-    """ A BaseImage instance encapsulating a Landsat-9 image.  Covers `region_*ha`.  """
+    """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`.  """
+    """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:
-    """ A BaseImage instance encapsulating a Landsat NDVI composite image.  Covers `region_*ha`.  """
+    """A BaseImage instance encapsulating a Landsat NDVI composite image.  Covers `region_*ha`."""
     return BaseImage.from_id(landsat_ndvi_image_id)
 
 
 @pytest.fixture(scope='session')
 def modis_nbar_base_image_unbnd(modis_nbar_image_id: str) -> BaseImage:
-    """ A BaseImage instance encapsulating an unbounded MODIS NBAR image in its native CRS.  Covers `region_*ha`.  """
+    """A BaseImage instance encapsulating an unbounded MODIS NBAR image in its native CRS.  Covers `region_*ha`."""
     return BaseImage(ee.Image(modis_nbar_image_id))
 
+
 @pytest.fixture(scope='session')
 def modis_nbar_base_image(modis_nbar_image_id: str, region_100ha: Dict) -> BaseImage:
-    """ A BaseImage instance encapsulating a MODIS NBAR image.  Covers `region_*ha`.  """
+    """A BaseImage instance encapsulating a MODIS NBAR image.  Covers `region_*ha`."""
     return BaseImage(ee.Image(modis_nbar_image_id).clip(region_100ha))
 
 
 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],
-    ]]
+    """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)))
@@ -130,7 +135,7 @@ def _test_export_profile(exp_profile, tgt_profile, transform_shape=False):
             assert exp_profile[key] == tgt_profile[key]
         assert exp_profile['transform'][:6] == pytest.approx(tgt_profile['transform'][:6], rel=1e-9)
     else:
-        assert exp_profile['transform'][:6] ==  pytest.approx(tgt_profile['transform'][:6], rel=0.05)
+        assert exp_profile['transform'][:6] == pytest.approx(tgt_profile['transform'][:6], rel=0.05)
 
     tgt_bounds = windows.bounds(
         windows.Window(0, 0, tgt_profile['width'], tgt_profile['height']), tgt_profile['transform']
@@ -145,16 +150,16 @@ def _test_export_profile(exp_profile, tgt_profile, transform_shape=False):
 
 
 def test_id_name(user_base_image: BaseImage, s2_sr_base_image: BaseImage):
-    """ Test `id` and `name` properties for different scenarios. """
+    """Test `id` and `name` properties for different scenarios."""
     assert user_base_image.id is None
     assert user_base_image.name is None
     # check that BaseImage.from_id() sets id without a getInfo
-    assert (s2_sr_base_image._id is not None)
+    assert s2_sr_base_image._id is not None
     assert s2_sr_base_image.name == s2_sr_base_image.id.replace('/', '-')
 
 
 def test_user_props(user_base_image: BaseImage):
-    """ Test non fixed projection image properties (other than id and has_fixed_projection). """
+    """Test non fixed projection image properties (other than id and has_fixed_projection)."""
     assert user_base_image.crs is None
     assert user_base_image.scale is None
     assert user_base_image.transform is None
@@ -167,7 +172,7 @@ def test_user_props(user_base_image: BaseImage):
 
 
 def test_fix_user_props(user_fix_base_image: BaseImage):
-    """ Test fixed projection image properties (other than id and has_fixed_projection). """
+    """Test fixed projection image properties (other than id and has_fixed_projection)."""
     assert user_fix_base_image.crs == 'EPSG:4326'
     assert user_fix_base_image.scale < 1
     assert user_fix_base_image.transform is not None
@@ -180,7 +185,7 @@ def test_fix_user_props(user_fix_base_image: BaseImage):
 
 
 def test_s2_props(s2_sr_base_image: BaseImage):
-    """ Test fixed projection S2 image properties (other than id and has_fixed_projection). """
+    """Test fixed projection S2 image properties (other than id and has_fixed_projection)."""
     min_band_info = s2_sr_base_image._ee_info['bands'][1]
     assert s2_sr_base_image.crs == min_band_info['crs']
     assert s2_sr_base_image.scale == min_band_info['crs_transform'][0]
@@ -198,7 +203,7 @@ def test_s2_props(s2_sr_base_image: BaseImage):
     'base_image', ['landsat_ndvi_base_image', 's2_sr_base_image', 'l9_base_image', 'modis_nbar_base_image']
 )
 def test_band_props(base_image: str, request: pytest.FixtureRequest):
-    """ Test `band_properties` completeness for generic/user/reflectance images. """
+    """Test `band_properties` completeness for generic/user/reflectance images."""
     base_image: BaseImage = request.getfixturevalue(base_image)
     assert base_image.band_properties is not None
     assert [bd['name'] for bd in base_image.band_properties] == [bd['id'] for bd in base_image._ee_info['bands']]
@@ -208,7 +213,7 @@ def test_band_props(base_image: str, request: pytest.FixtureRequest):
 
 
 def test_has_fixed_projection(user_base_image: BaseImage, user_fix_base_image: BaseImage, s2_sr_base_image):
-    """ Test the `has_fixed_projection` property.  """
+    """Test the `has_fixed_projection` property."""
     assert not user_base_image.has_fixed_projection
     assert user_fix_base_image.has_fixed_projection
     assert s2_sr_base_image.has_fixed_projection
@@ -230,7 +235,7 @@ def test_has_fixed_projection(user_base_image: BaseImage, user_fix_base_image: B
 )
 # yapf: enable
 def test_min_dtype(ee_data_type_list: List, exp_dtype: str):
-    """ Test BasicImage.__get_min_dtype() with emulated EE info dicts.  """
+    """Test BasicImage.__get_min_dtype() with emulated EE info dicts."""
     ee_info = dict(bands=[])
     for ee_data_type in ee_data_type_list:
         ee_info['bands'].append(dict(data_type=ee_data_type))
@@ -238,19 +243,20 @@ def test_min_dtype(ee_data_type_list: List, exp_dtype: str):
 
 
 def test_convert_dtype_error():
-    """ Test BaseImage.test_convert_dtype() raises an error with incorrect dtype. """
+    """Test BaseImage.test_convert_dtype() raises an error with incorrect dtype."""
     with pytest.raises(TypeError):
         BaseImage._convert_dtype(ee.Image(1), dtype='unknown')
 
 
 @pytest.mark.parametrize('size, exp_str', [(1024, '1.02 KB'), (234.56e6, '234.56 MB'), (1e9, '1.00 GB')])
 def test_str_format_size(size: int, exp_str: str):
-    """ Test formatting of byte sizes as human readable strings. """
+    """Test formatting of byte sizes as human readable strings."""
     assert BaseImage._str_format_size(size) == exp_str
 
 
 @pytest.mark.parametrize(
-    'params', [
+    'params',
+    [
         dict(),
         dict(crs='EPSG:4326'),
         dict(crs='EPSG:4326', region='region_25ha'),
@@ -259,7 +265,7 @@ def test_str_format_size(size: int, exp_str: str):
         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):
     """
@@ -274,14 +280,15 @@ def test_prepare_no_fixed_projection(user_base_image: BaseImage, params: Dict, r
 
 
 @pytest.mark.parametrize(
-    'base_image, params', [
+    'base_image, params',
+    [
         ('user_fix_base_image', dict()),
         ('modis_nbar_base_image_unbnd', dict(crs='EPSG:4326')),
         ('modis_nbar_base_image_unbnd', dict(crs='EPSG:4326', scale=100)),
         ('modis_nbar_base_image_unbnd', dict(crs='EPSG:4326', crs_transform=Affine.identity())),
         ('modis_nbar_base_image_unbnd', dict(crs='EPSG:4326', shape=(100, 100))),
         ('modis_nbar_base_image_unbnd', dict(crs_transform=Affine.identity(), shape=(100, 100))),
-    ]
+    ],
 )  # yapf: disable
 def test_prepare_unbounded(base_image: BaseImage, params: Dict, request: pytest.FixtureRequest):
     """
@@ -295,7 +302,7 @@ def test_prepare_unbounded(base_image: BaseImage, params: Dict, request: pytest.
 
 
 def test_prepare_exceptions(user_base_image: BaseImage, user_fix_base_image: BaseImage, region_25ha: Dict):
-    """ Test remaining BaseImage._prepare_for_export() error cases. """
+    """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)
@@ -307,10 +314,11 @@ def test_prepare_exceptions(user_base_image: BaseImage, user_fix_base_image: Bas
 
 
 @pytest.mark.parametrize(
-    'base_image', ['s2_sr_base_image', 'l9_base_image', 'modis_nbar_base_image'],
+    'base_image',
+    ['s2_sr_base_image', 'l9_base_image', 'modis_nbar_base_image'],
 )
 def test_prepare_defaults(base_image: str, request: pytest.FixtureRequest):
-    """ Test BaseImage._prepare_for_export() with no (i.e. default) arguments with bounded images. """
+    """Test BaseImage._prepare_for_export() with no (i.e. default) arguments with bounded images."""
     base_image: BaseImage = request.getfixturevalue(base_image)
     exp_image = base_image._prepare_for_export()
 
@@ -320,7 +328,8 @@ def test_prepare_defaults(base_image: str, request: pytest.FixtureRequest):
 
 
 @pytest.mark.parametrize(
-    'base_image, param_image', [
+    'base_image, param_image',
+    [
         ('s2_sr_base_image', 's2_sr_base_image'),
         ('l9_base_image', 's2_sr_base_image'),
         ('modis_nbar_base_image', 's2_sr_base_image'),
@@ -336,10 +345,10 @@ def test_prepare_defaults(base_image: str, request: pytest.FixtureRequest):
         ('modis_nbar_base_image', 'modis_nbar_base_image'),
         ('user_base_image', 'modis_nbar_base_image'),
         ('user_fix_base_image', 'modis_nbar_base_image'),
-    ]  # yapf: disable
+    ],  # yapf: disable
 )
 def test_prepare_transform_shape(base_image: str, param_image: str, request: pytest.FixtureRequest):
-    """ Test BaseImage._prepare_for_export() with crs_transform and shape parameters. """
+    """Test BaseImage._prepare_for_export() with crs_transform and shape parameters."""
     base_image: BaseImage = request.getfixturevalue(base_image)
     param_image: BaseImage = request.getfixturevalue(param_image)
 
@@ -356,7 +365,8 @@ def test_prepare_transform_shape(base_image: str, param_image: str, request: pyt
 
 
 @pytest.mark.parametrize(
-    'base_image, param_image', [
+    'base_image, param_image',
+    [
         ('s2_sr_base_image', 's2_sr_base_image'),
         ('l9_base_image', 's2_sr_base_image'),
         ('modis_nbar_base_image', 's2_sr_base_image'),
@@ -372,13 +382,13 @@ def test_prepare_transform_shape(base_image: str, param_image: str, request: pyt
         ('modis_nbar_base_image', 'modis_nbar_base_image'),
         ('user_base_image', 'modis_nbar_base_image'),
         ('user_fix_base_image', 'modis_nbar_base_image'),
-    ]
+    ],
 )  # yapf: disable
 def test_prepare_region_scale(base_image: str, param_image: str, region_25ha: dict, request: pytest.FixtureRequest):
-    """ Test BaseImage._prepare_for_export() with region and scale parameters. """
+    """Test BaseImage._prepare_for_export() with region and scale parameters."""
     base_image: BaseImage = request.getfixturevalue(base_image)
     _param_image: BaseImage = request.getfixturevalue(param_image)
-    param_image = copy.deepcopy(_param_image)   # avoid changing session fixture
+    param_image = copy.deepcopy(_param_image)  # avoid changing session fixture
     # clip param_image, so that param_image.transform can be tested against below
     param_image.ee_image = param_image.ee_image.clip(region_25ha)
 
@@ -402,19 +412,22 @@ def test_prepare_region_scale(base_image: str, param_image: str, region_25ha: di
     exp_crs = CRS.from_string(utils.rio_crs(exp_image.crs))
     tgt_bounds = warp.transform_bounds(tgt_footprint_crs, exp_crs, *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])
+        (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])
     )
 
 
 @pytest.mark.parametrize(
-    'base_image, bands', [
+    'base_image, bands',
+    [
         ('s2_sr_base_image', ['B1', 'B5']),
         ('l9_base_image', ['SR_B4', 'SR_B3', 'SR_B2']),
-    ]
+    ],
 )  # yapf: disable
 def test_prepare_bands(base_image: str, bands: List[str], region_25ha: dict, request: pytest.FixtureRequest):
-    """ Test BaseImage._prepare_for_export() with bands parameter. """
+    """Test BaseImage._prepare_for_export() with bands parameter."""
     base_image: BaseImage = request.getfixturevalue(base_image)
     param_image = BaseImage(base_image.ee_image.select(bands))
 
@@ -426,18 +439,17 @@ def test_prepare_bands(base_image: str, bands: List[str], region_25ha: dict, req
 
 
 def test_prepare_bands_error(s2_sr_base_image):
-    """ Test BaseImage._prepare_for_export() raises an error with incorrect bands. """
+    """Test BaseImage._prepare_for_export() raises an error with incorrect bands."""
     with pytest.raises(ValueError):
         s2_sr_base_image._prepare_for_export(bands=['unknown'])
 
 
 @pytest.mark.parametrize(
-    'base_image', ['s2_sr_base_image', 'l9_base_image', 'modis_nbar_base_image'],
+    'base_image',
+    ['s2_sr_base_image', 'l9_base_image', 'modis_nbar_base_image'],
 )
-def test_prepare_for_download(
-    base_image: str, request: pytest.FixtureRequest
-):
-    """ Test BaseImage._prepare_for_download() sets rasterio profile as expected.  """
+def test_prepare_for_download(base_image: str, request: pytest.FixtureRequest):
+    """Test BaseImage._prepare_for_download() sets rasterio profile as expected."""
     base_image: BaseImage = request.getfixturevalue(base_image)
     exp_image, exp_profile = base_image._prepare_for_download()
     tgt_profile = base_image.profile
@@ -446,13 +458,20 @@ def test_prepare_for_download(
 
 
 @pytest.mark.parametrize(
-    'dtype, exp_nodata', [
-        ('uint8', 0), ('int8', -2**7), ('uint16', 0), ('int16', -2**15), ('uint32', 0), ('int32', -2**31),
-        ('float32', float('nan')), ('float64', float('nan')),
-    ]
+    'dtype, exp_nodata',
+    [
+        ('uint8', 0),
+        ('int8', -(2**7)),
+        ('uint16', 0),
+        ('int16', -(2**15)),
+        ('uint32', 0),
+        ('int32', -(2**31)),
+        ('float32', float('nan')),
+        ('float64', float('nan')),
+    ],
 )
 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.  """
+    """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']):
@@ -462,14 +481,15 @@ def test_prepare_nodata(user_fix_base_image: BaseImage, region_25ha: Dict, dtype
 
 
 @pytest.mark.parametrize(
-    'src_image, dtype', [
+    'src_image, dtype',
+    [
         ('s2_sr_base_image', 'float32'),
         ('l9_base_image', None),
         ('modis_nbar_base_image', None),
-    ]
+    ],
 )  # yapf: disable
 def test_scale_offset(src_image: str, dtype: str, region_100ha: Dict, request: pytest.FixtureRequest):
-    """ Test BaseImage._prepare_for_export(scale_offset=True) gives expected properties and reflectance ranges. """
+    """Test BaseImage._prepare_for_export(scale_offset=True) gives expected properties and reflectance ranges."""
 
     src_image: BaseImage = request.getfixturevalue(src_image)
     exp_image = src_image._prepare_for_export(scale_offset=True)
@@ -479,7 +499,7 @@ def test_scale_offset(src_image: str, dtype: str, region_100ha: Dict, request: p
     assert exp_image.dtype == dtype or 'float64'
 
     def get_min_max_refl(base_image: BaseImage) -> Dict:
-        """ Get the min & max of each reflectance band of base_image.  """
+        """Get the min & max of each reflectance band of base_image."""
         band_props = base_image.band_properties
         refl_bands = [
             bp['name'] for bp in band_props if ('center_wavelength' in bp) and (bp['center_wavelength'] < 1)
@@ -487,7 +507,7 @@ def get_min_max_refl(base_image: BaseImage) -> Dict:
         ee_image = base_image.ee_image.select(refl_bands)
         min_max_dict = ee_image.reduceRegion(
             reducer=ee.Reducer.minMax(), geometry=region_100ha, bestEffort=True
-        ).getInfo() # yapf: disable
+        ).getInfo()  # yapf: disable
         min_dict = {k: v for k, v in min_max_dict.items() if 'min' in k}
         max_dict = {k: v for k, v in min_max_dict.items() if 'max' in k}
         return min_dict, max_dict
@@ -499,7 +519,7 @@ def get_min_max_refl(base_image: BaseImage) -> Dict:
 
 
 def test_tile_shape():
-    """ Test BaseImage._get_tile_shape() satisfies the tile size limit for different image shapes. """
+    """Test BaseImage._get_tile_shape() satisfies the tile size limit for different image shapes."""
     max_tile_dim = 10000
 
     for max_tile_size in range(4, 32, 4):
@@ -517,14 +537,15 @@ def test_tile_shape():
 
 
 @pytest.mark.parametrize(
-    'image_shape, tile_shape, image_transform', [
+    'image_shape, tile_shape, image_transform',
+    [
         ((1000, 500), (101, 101), Affine.identity()),
         ((1000, 100), (101, 101), Affine.scale(1.23)),
         ((1000, 102), (101, 101), Affine.scale(1.23) * Affine.translation(12, 34)),
-    ]
+    ],
 )  # yapf: disable
 def test_tiles(image_shape: Tuple, tile_shape: Tuple, image_transform: Affine):
-    """ Test continuity and coverage of tiles. """
+    """Test continuity and coverage of tiles."""
     exp_image = BaseImageLike(shape=image_shape, transform=image_transform)
     tiles = [tile for tile in exp_image._tiles(tile_shape=tile_shape)]
 
@@ -550,14 +571,12 @@ def test_tiles(image_shape: Tuple, tile_shape: Tuple, image_transform: Affine):
 
 
 def test_download_transform_shape(user_base_image: str, tmp_path: pathlib.Path, request: pytest.FixtureRequest):
-    """ Test download file properties and pixel data with crs_transform and shape arguments. """
+    """Test download file properties and pixel data with crs_transform and shape arguments."""
     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']
-    )
+    tgt_bounds = windows.bounds(windows.Window(0, 0, *shape[::-1]), tgt_prof['transform'])
     download_params = dict(
         crs=tgt_prof['crs'], crs_transform=tgt_prof['transform'], shape=shape, dtype=tgt_prof['dtype']
     )
@@ -572,14 +591,12 @@ def test_download_transform_shape(user_base_image: str, tmp_path: pathlib.Path,
 
 
 def test_download_region_scale(user_base_image: str, tmp_path: pathlib.Path, request: pytest.FixtureRequest):
-    """ Test download file properties and pixel data with region and scale arguments. """
+    """Test download file properties and pixel data with region and scale arguments."""
     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']
-    )
+    tgt_bounds = windows.bounds(windows.Window(0, 0, *shape[::-1]), tgt_prof['transform'])
     region = bounds_polygon(*tgt_bounds, crs=tgt_prof['crs'])
     download_params = dict(crs=tgt_prof['crs'], region=region, scale=tgt_prof['transform'][0], dtype=tgt_prof['dtype'])
     filename = tmp_path.joinpath('test.tif')
@@ -593,7 +610,7 @@ def test_download_region_scale(user_base_image: str, tmp_path: pathlib.Path, req
 
 
 def test_overviews(user_base_image: BaseImage, region_25ha: Dict, tmp_path: pathlib.Path):
-    """ Test overviews get built on download. """
+    """Test overviews get built on download."""
     filename = tmp_path.joinpath('test_user_download.tif')
     user_base_image.download(filename, region=region_25ha, crs='EPSG:3857', scale=1)
     assert filename.exists()
@@ -604,7 +621,7 @@ def test_overviews(user_base_image: BaseImage, region_25ha: Dict, tmp_path: path
 
 
 def test_metadata(s2_sr_base_image: BaseImage, region_25ha: Dict, tmp_path: pathlib.Path):
-    """ Test metadata is written to a downloaded file. """
+    """Test metadata is written to a downloaded file."""
     filename = tmp_path.joinpath('test_s2_band_subset_download.tif')
     s2_sr_base_image.download(filename, region=region_25ha, crs='EPSG:3857', scale=60, bands=['B9'])
     assert filename.exists()
@@ -619,18 +636,18 @@ def test_metadata(s2_sr_base_image: BaseImage, region_25ha: Dict, tmp_path: path
 
 @pytest.mark.parametrize('type', [ExportType.drive, ExportType.asset, ExportType.cloud])  # yapf: disable
 def test_start_export(type, user_fix_base_image: BaseImage, region_25ha: Dict):
-    """ Test start of a small export. """
+    """Test start of a small export."""
     # Note that this export should start successfully, but will ultimately fail for the asset and cloud options.  For
     # the asset option, there will be an overwrite issue.  For cloud storage, there is no 'geedim' bucket.
     task = user_fix_base_image.export(
         'test_export', type=type, folder='geedim', scale=30, region=region_25ha, wait=False
     )
     assert task.active()
-    assert task.status()['state'] == 'READY'
+    assert task.status()['state'] == 'READY' or task.status()['state'] == 'RUNNING'
 
 
 def __test_export_asset(user_fix_base_image: BaseImage, region_25ha: Dict):
-    """ Test start of a small export to EE asset. """
+    """Test start of a small export to EE asset."""
     # TODO: consider removing this slow test in favour of the equivalent integration test
     filename = f'test_export_{np.random.randint(1<<31)}'
     folder = 'geedim'
@@ -655,7 +672,7 @@ def __test_export_asset(user_fix_base_image: BaseImage, region_25ha: Dict):
 
 
 def test_download_bigtiff(s2_sr_base_image: BaseImage, tmp_path: pathlib.Path):
-    """ Test that BIGTIFF gets set in the profile of images larger than 4GB. """
+    """Test that BIGTIFF gets set in the profile of images larger than 4GB."""
     exp_image, profile = s2_sr_base_image._prepare_for_download()
     assert exp_image.size >= 4e9
     assert 'bigtiff' in profile
@@ -663,23 +680,24 @@ def test_download_bigtiff(s2_sr_base_image: BaseImage, tmp_path: pathlib.Path):
 
 
 def test_prepare_ee_geom(l9_base_image: BaseImage, tmp_path: pathlib.Path):
-    """ Test that _prepare_for_export works with an ee.Geometry region at native crs and scale (Issue #6). """
+    """Test that _prepare_for_export works with an ee.Geometry region at native crs and scale (Issue #6)."""
     region = l9_base_image.ee_image.geometry()
     exp_image = l9_base_image._prepare_for_export(region=region)
     assert exp_image.scale == l9_base_image.scale
 
 
 @pytest.mark.parametrize(
-    'base_image, exp_value', [
+    'base_image, exp_value',
+    [
         ('s2_sr_base_image', True),
         ('l9_base_image', True),
         ('modis_nbar_base_image_unbnd', False),
         ('user_base_image', False),
         ('user_fix_base_image', False),
-    ]
+    ],
 )  # yapf: disable
 def test_bounded(base_image: str, exp_value: bool, request: pytest.FixtureRequest):
-    """ Test BaseImage.bounded has correct value for different images. """
+    """Test BaseImage.bounded has correct value for different images."""
     base_image: BaseImage = request.getfixturevalue(base_image)
     assert base_image.bounded == exp_value
 
diff --git a/tests/test_tile.py b/tests/test_tile.py
index dcece64..7739c96 100644
--- a/tests/test_tile.py
+++ b/tests/test_tile.py
@@ -13,23 +13,32 @@
     See the License for the specific language governing permissions and
     limitations under the License.
 """
+
+import io
+import json
+import logging
+import zipfile
 from collections import namedtuple
 
 import ee
 import numpy as np
 import pytest
-from geedim.tile import Tile
-from geedim.utils import retry_session
+import rasterio as rio
+import requests
 from rasterio import Affine
 from rasterio.windows import Window
 from tqdm.auto import tqdm
 
+from geedim.errors import TileError
+from geedim.tile import Tile
+from geedim.utils import retry_session
+
 BaseImageLike = namedtuple('BaseImageLike', ['ee_image', 'crs', 'transform', 'shape', 'count', 'dtype'])
 
 
 @pytest.fixture(scope='module')
-def base_image_like(region_25ha):
-    """ Create a synthetic image object to emulate BaseImage. """
+def mock_base_image(region_25ha: dict) -> BaseImageLike:
+    """A BaseImage mock containing a synthetic ee.Image."""
     ee_image = ee.Image([1, 2, 3]).reproject(crs='EPSG:4326', scale=30).clip(region_25ha)
     ee_info = ee_image.getInfo()
     band_info = ee_info['bands'][0]
@@ -37,28 +46,155 @@ def base_image_like(region_25ha):
     return BaseImageLike(ee_image, 'EPSG:3857', transform, tuple(band_info['dimensions'][::-1]), 3, 'uint8')
 
 
-def test_create(base_image_like):
-    """ Test creation of a Tile object that refers to the whole of `base_image_like`. """
-    window = Window(0, 0, *base_image_like.shape[::-1])
-    tile = Tile(base_image_like, window)
+@pytest.fixture(scope='module')
+def synth_tile(mock_base_image: BaseImageLike) -> Tile:
+    """A tile representing the whole of ``mock_base_image``."""
+    window = Window(0, 0, *mock_base_image.shape[::-1])
+    return Tile(mock_base_image, window)
+
+
+@pytest.fixture(scope='function')
+def mock_ee_image(monkeypatch: pytest.MonkeyPatch) -> None:
+    """Patch ee.Image.getDownloadURL() to return None."""
+
+    def getDownloadURL(*args, **kwargs):
+        return None
+
+    monkeypatch.setattr(ee.Image, 'getDownloadURL', getDownloadURL)
+
+
+@pytest.fixture(scope='module')
+def zipped_gtiff_bytes(mock_base_image: BaseImageLike) -> bytes:
+    """Zipped GeoTIFF bytes for ``mock_base_image``."""
+    zip_buffer = io.BytesIO()
+    array = np.ones((mock_base_image.count, *mock_base_image.shape)) * np.array([1, 2, 3]).reshape(-1, 1, 1)
+
+    with rio.MemoryFile() as mem_file:
+        with mem_file.open(
+            **rio.default_gtiff_profile,
+            width=mock_base_image.shape[1],
+            height=mock_base_image.shape[0],
+            count=mock_base_image.count
+        ) as ds:
+            ds.write(array)
+        with zipfile.ZipFile(zip_buffer, 'a', zipfile.ZIP_DEFLATED, False) as zf:
+            zf.writestr('test.tif', mem_file.read())
+
+    return zip_buffer.getvalue()
+
+
+def test_create(mock_base_image: BaseImageLike):
+    """Test creation of a Tile object."""
+    window = Window(0, 0, *mock_base_image.shape[::-1])
+    tile = Tile(mock_base_image, window)
     assert tile.window == window
-    assert tile._transform == base_image_like.transform
-    assert tile._shape == base_image_like.shape
+    assert tile._transform == mock_base_image.transform
+    assert tile._shape == mock_base_image.shape
 
 
 @pytest.mark.parametrize('session', [None, retry_session()])
-def test_download(base_image_like, session):
-    """ Test downloading the synthetic image tile.  """
-    window = Window(0, 0, *base_image_like.shape[::-1])
-    tile = Tile(base_image_like, window)
-    dtype_size = np.dtype(tile._exp_image.dtype).itemsize
-    raw_download_size = tile._shape[0] * tile._shape[1] * tile._exp_image.count * dtype_size
+def test_download(synth_tile: Tile, session):
+    """Test downloading the synthetic image tile."""
+    dtype_size = np.dtype(synth_tile._exp_image.dtype).itemsize
+    raw_download_size = synth_tile._shape[0] * synth_tile._shape[1] * synth_tile._exp_image.count * dtype_size
     bar = tqdm(total=float(raw_download_size))
-    array = tile.download(session=session, bar=bar)
+    array = synth_tile.download(session=session, bar=bar)
 
     assert array is not None
-    assert array.shape == (base_image_like.count, *base_image_like.shape)
-    assert array.dtype == np.dtype(base_image_like.dtype)
-    for i in range(3):
+    assert array.shape == (synth_tile._exp_image.count, *synth_tile._exp_image.shape)
+    assert array.dtype == np.dtype(synth_tile._exp_image.dtype)
+    for i in range(array.shape[0]):
         assert np.all(array[i] == i + 1)
     assert bar.n == pytest.approx(raw_download_size, rel=0.01)
+
+
+def test_mem_limit_error(synth_tile: Tile, mock_ee_image: None):
+    """Test downloading raises the 'user memory limit exceeded' error with a mock response."""
+    # patch session.get() to return a mock response with EE memory limit error
+    session = retry_session()
+    msg = 'User memory limit exceeded.'
+
+    def get(url, **kwargs):
+        response = requests.Response()
+        response.status_code = 400
+        response.headers = {'content-length': '1'}
+        response._content = json.dumps({'error': {'message': msg}}).encode()
+        return response
+
+    session.get = get
+
+    # test memory limit error is raised on download
+    with pytest.raises(TileError) as ex:
+        synth_tile.download(session=session)
+    assert msg in str(ex.value)
+
+
+def test_retry(synth_tile: Tile, mock_ee_image: None, zipped_gtiff_bytes: bytes, caplog: pytest.LogCaptureFixture):
+    """Test downloading retries invalid tiles until it succeeds."""
+    # create progress bar
+    dtype_size = np.dtype(synth_tile._exp_image.dtype).itemsize
+    raw_download_size = synth_tile._shape[0] * synth_tile._shape[1] * synth_tile._exp_image.count * dtype_size
+    bar = tqdm(total=float(raw_download_size))
+
+    # create mock invalid responses for each retry
+    responses = []
+    for _ in range(5):
+        response = requests.Response()
+        response.status_code = 200
+        response.headers = {'content-length': str(len(zipped_gtiff_bytes))}
+        response.raw = io.BytesIO(b'error')
+        responses.append(response)
+
+    # make the last response valid
+    responses[-1].raw = io.BytesIO(zipped_gtiff_bytes)
+
+    # patch session.get() to pop and return a mocked response from the list
+    session = retry_session()
+
+    def get(url, **kwargs):
+        return responses.pop(0)
+
+    session.get = get
+
+    # test the tile is downloaded correctly, after retries
+    with caplog.at_level(logging.WARNING):
+        array = synth_tile.download(session=session, bar=bar, backoff_factor=0)
+
+    assert array.shape == (synth_tile._exp_image.count, *synth_tile._exp_image.shape)
+    assert array.dtype == np.dtype(synth_tile._exp_image.dtype)
+    for i in range(array.shape[0]):
+        assert np.all(array[i] == i + 1)
+
+    # test progress bar is adjusted for retries
+    assert bar.n == pytest.approx(raw_download_size, rel=0.01)
+
+    # test retry logs
+    assert 'retry' in caplog.text and 'zip' in caplog.text
+
+
+def test_retry_error(synth_tile: Tile, mock_ee_image: None, zipped_gtiff_bytes: bytes):
+    """Test downloading raises an error when the maximum retries are reached."""
+    # create progress bar
+    dtype_size = np.dtype(synth_tile._exp_image.dtype).itemsize
+    raw_download_size = synth_tile._shape[0] * synth_tile._shape[1] * synth_tile._exp_image.count * dtype_size
+    bar = tqdm(total=float(raw_download_size))
+
+    # patch session.get() to return a mock response with invalid bytes
+    session = retry_session()
+
+    def get(url, **kwargs):
+        response = requests.Response()
+        response.status_code = 200
+        response.headers = {'content-length': '10'}
+        response.raw = io.BytesIO(b'error')
+        return response
+
+    session.get = get
+
+    # test max retries error is raised on download
+    with pytest.raises(TileError) as ex:
+        synth_tile.download(session=session, bar=bar, backoff_factor=0)
+    assert 'maximum retries' in str(ex.value)
+
+    # test progress bar is adjusted for retries
+    assert bar.n == pytest.approx(0)