tifffile

1
# tifffile/examples/issues125.py
2

3
"""Create a Fsspec ReferenceFileSystem for a sequence of TIFF files on S3
4

5
This Python script uses the Tifffile and Fsspec libraries to create a
6
multiscale ReferenceFileSystem JSON file for a sequence of cloud optimized
7
GeoTIFF (COG) files stored on S3. The tiles of the COG files are used as
8
chunks. No additional Numcodecs codec needs to be registered since the COG
9
files use Zlib compression. A Xarray dataset is created from the
10
ReferenceFileSystem file and a subset of the dataset is plotted.
11

12
See https://github.com/cgohlke/tifffile/issues/125
13

14
"""
15

16
import fsspec
17
import tifffile
18
import xarray
19
from matplotlib import pyplot
20

21
# get a list of cloud optimized GeoTIFF files stored on S3
22
remote_options = {
23
    'anon': True,
24
    'client_kwargs': {'endpoint_url': 'https://mghp.osn.xsede.org'},
25
}
26
fs = fsspec.filesystem('s3', **remote_options)
27
files = [f's3://{f}' for f in fs.ls('/rsignellbucket1/lcmap/cog')]
28

29
# write the ReferenceFileSystem of each file to a JSON file
30
with open('issue125.json', 'w', encoding='utf-8', newline='\n') as jsonfile:
31
    for i, filename in enumerate(tifffile.natural_sorted(files)):
32
        url, name = filename.rsplit('/', 1)
33
        with fs.open(filename) as fh:
34
            with tifffile.TiffFile(fh, name=name) as tif:
35
                print(tif.geotiff_metadata)
36
                with tif.series[0].aszarr() as store:
37
                    store.write_fsspec(
38
                        jsonfile,
39
                        url=url,
40
                        # using an experimental API:
41
                        _shape=[len(files)],  # shape of file sequence
42
                        _axes='T',  # axes of file sequence
43
                        _index=[i],  # index of this file in sequence
44
                        _append=i != 0,  # if True, only write index keys+value
45
                        _close=i == len(files) - 1,  # if True, no more appends
46
                        # groupname='0',  # required for non-pyramidal series
47
                    )
48

49
# create a fsspec mapper instance from the ReferenceFileSystem file
50
mapper = fsspec.get_mapper(
51
    'reference://',
52
    fo='issue125.json',
53
    target_protocol='file',
54
    remote_protocol='s3',
55
    remote_options=remote_options,
56
)
57

58
# create a xarray dataset from the mapper
59
dataset = xarray.open_zarr(mapper, consolidated=False)
60
print(dataset)
61

62
# plot a slice of the 5th pyramidal level of the dataset
63
xarray.plot.imshow(dataset['5'][0, 32:-32, 32:-32])
64
pyplot.show()
65