Introduction

This notebook reads in the red, blue and green bands from a landsat 8 image and uses rasterio to write out new 3-channel color tiff, jpeg and png files

In [1]:
from pathlib import Path
from a301.utils.data_read import download
import numpy as np
import pprint

Get the vancouver image

In [2]:
import rasterio
import a301
filenames=["LC08_L1TP_047026_20150614_20180131_01_T1_B2.TIF",
    "LC08_L1TP_047026_20150614_20180131_01_T1_B3.TIF",
    "LC08_L1TP_047026_20150614_20180131_01_T1_B4.TIF",
    "LC08_L1TP_047026_20150614_20180131_01_T1_MTL.txt"]
dest_folder=a301.data_dir / Path("landsat8/vancouver")
In [3]:
for the_file in filenames:
    landsat_tif = Path('landsat_scenes/l8_vancouver') / Path(the_file)
    download(str(landsat_tif),dest_folder=dest_folder)
band2=list(dest_folder.glob("*_B2.TIF"))[0]
band3=list(dest_folder.glob("*_B3.TIF"))[0]
band4=list(dest_folder.glob("*_B4.TIF"))[0]
trying https://clouds.eos.ubc.ca/~phil/courses/atsc301/downloads/landsat_scenes/l8_vancouver/LC08_L1TP_047026_20150614_20180131_01_T1_B2.TIF
writing to: C:\Users\phil\repos\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_B2.TIF
writing temporary file C:\Users\phil\repos\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_B2.TIF_tmp
downloaded LC08_L1TP_047026_20150614_20180131_01_T1_B2.TIF
size = 70873586
trying https://clouds.eos.ubc.ca/~phil/courses/atsc301/downloads/landsat_scenes/l8_vancouver/LC08_L1TP_047026_20150614_20180131_01_T1_B3.TIF
writing to: C:\Users\phil\repos\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_B3.TIF
writing temporary file C:\Users\phil\repos\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_B3.TIF_tmp
downloaded LC08_L1TP_047026_20150614_20180131_01_T1_B3.TIF
size = 74130783
trying https://clouds.eos.ubc.ca/~phil/courses/atsc301/downloads/landsat_scenes/l8_vancouver/LC08_L1TP_047026_20150614_20180131_01_T1_B4.TIF
writing to: C:\Users\phil\repos\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_B4.TIF
writing temporary file C:\Users\phil\repos\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_B4.TIF_tmp
downloaded LC08_L1TP_047026_20150614_20180131_01_T1_B4.TIF
size = 75252265
trying https://clouds.eos.ubc.ca/~phil/courses/atsc301/downloads/landsat_scenes/l8_vancouver/LC08_L1TP_047026_20150614_20180131_01_T1_MTL.txt
writing to: C:\Users\phil\repos\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_MTL.txt
writing temporary file C:\Users\phil\repos\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_MTL.txt_tmp
downloaded LC08_L1TP_047026_20150614_20180131_01_T1_MTL.txt
size = 8702
Warning -- your file is tiny (smaller than 10 Kbyte)
Did something go wrong?

Read in the bands 2 (blue), 3 (green) and 4 (red)

Note that rasterio is a pretty complicated object with a lot of functionality.

The full documentation is at https://rasterio.readthedocs.io/en/latest/

Save the image profile, coordinate reference system (crs) and affine transform for inspection -- these will be the same for every band in the image

In [4]:
with rasterio.open(band2) as b2_raster:
    b2_data = b2_raster.read(1)
    transform=b2_raster.transform
    crs=b2_raster.crs
    profile=b2_raster.profile
with rasterio.open(band3) as b3_raster:
    b3_data = b3_raster.read(1)
with rasterio.open(band4) as b4_raster:
    b4_data = b4_raster.read(1)

The profile for this image

In [5]:
print(f"\nprofile: {pprint.pformat(profile)}\n")
profile: {'blockxsize': 256,
 'blockysize': 256,
 'compress': 'lzw',
 'count': 1,
 'crs': CRS({'init': 'epsg:32610'}),
 'driver': 'GTiff',
 'dtype': 'uint16',
 'height': 7961,
 'interleave': 'band',
 'nodata': None,
 'tiled': True,
 'transform': Affine(30.0, 0.0, 397185.0,
       0.0, -30.0, 5531415.0),
 'width': 7851}

Create new 3-dimensional array and copy the bands

Make an empty 3-d array to hold the three channels. We are going to scale each band to the range 0-255 so specify 8 bit np.uint8 words

In [6]:
channels=np.empty([3,b2_data.shape[0],b2_data.shape[1]],dtype=np.uint8)

Stretch and write out the combined image

Do a histogram stretch on each band using skimage.exposure and save it into the 3-d array.

In [7]:
from skimage import  img_as_ubyte
from skimage import exposure
for index,image in enumerate([b4_data,b3_data,b2_data]):
    stretched=exposure.equalize_hist(image)    
    channels[index,:,:] = img_as_ubyte(stretched)
C:\Users\phil\ma36\envs\fiona\lib\site-packages\skimage\util\dtype.py:141: UserWarning: Possible precision loss when converting from float64 to uint8
  .format(dtypeobj_in, dtypeobj_out))

Write out the geotiff file

In [8]:
tif_filename = dest_folder / Path('vancouver_432.tiff')
num_chans, height, width = channels.shape
with rasterio.open(tif_filename,'w',driver='GTiff',
                   height=height,width=width,
                   count=num_chans,dtype=channels.dtype,
                   crs=crs,transform=transform, nodata=0.0) as dst:
        dst.write(channels)
        keys=['4','3','2']
        for index,chan_name in enumerate(keys):
            dst.update_tags(index+1,name=chan_name)

Read the geotiff file and convert to png and jpeg

In [9]:
with rasterio.open(tif_filename) as infile:
    print(f"\nnew profile: {pprint.pformat(infile.profile)}\n")
    profile=infile.profile
    #
    # change the driver name from GTiff to PNG
    #
    profile['driver']='PNG'
    #
    # pathlib makes it easy to add a new suffix to a
    # filename
    #    
    png_filename=tif_filename.with_suffix('.png')
    raster=infile.read()
    with rasterio.open(png_filename, 'w', **profile) as dst:
        dst.write(raster)
    #
    # now do jpeg
    #
    profile['driver']='JPEG'
    jpeg_filename=tif_filename.with_suffix('.jpeg')
    with rasterio.open(jpeg_filename, 'w', **profile) as dst:
        dst.write(raster)
new profile: {'count': 3,
 'crs': CRS({'init': 'epsg:32610'}),
 'driver': 'GTiff',
 'dtype': 'uint8',
 'height': 7961,
 'interleave': 'pixel',
 'nodata': 0.0,
 'tiled': False,
 'transform': Affine(30.0, 0.0, 397185.0,
       0.0, -30.0, 5531415.0),
 'width': 7851}

Display the jpeg file

In [10]:
from IPython.display import Image
Image(str(jpeg_filename))
Out[10]: