<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Introduction2&nbsp;&nbsp;Get bands 3, 4  (green, red)3&nbsp;&nbsp;Rename the band 5 tiff file created by the image_zoom notebook4&nbsp;&nbsp;Open the small file and get the transform5&nbsp;&nbsp;Use the transform to make a Window with the correct ul corner6&nbsp;&nbsp;Read only the window pixels from the band 3, 4 files7&nbsp;&nbsp;check to make sure the read worked8&nbsp;&nbsp;Compare the green and red bands9&nbsp;&nbsp;Write out a new composite  IntroductionWe've got a new tiff file for band 5 called small_file.tiff that was written out by the image_zoom notebook.  We can use the corners of that image to clip other bands and saveus the trouble of working with the large original files.Below I use the [rasterio windows module](https://rasterio.readthedocs.io/en/latest/topics/windowed-rw.html) to read in only the part of the image that matches small_file.tiff<jupyter_code>import rasterio
import a301
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import Normalize
from a301.landsat.landsat_metadata import landsat_metadata
import cartopy
from rasterio.windows import Window
from pyproj import transform as proj_transform
from pyproj import Proj
from a301.landsat.toa_reflectance import calc_reflc_8
import pprint
from a301.utils.data_read import download
from pathlib import Path
from affine import Affine
from IPython.display import Image<jupyter_output><empty_output><jupyter_text>Get bands 3, 4  (green, red)<jupyter_code>filenames=["LC08_L1TP_047026_20150614_20180131_01_T1_B4.TIF",
    "LC08_L1TP_047026_20150614_20180131_01_T1_B3.TIF",
    "LC08_L1TP_047026_20150614_20180131_01_T1_MTL.txt"]
dest_folder=a301.data_dir / Path("landsat8/vancouver")
dest_folder.mkdir(parents=True, exist_ok=True)
for the_file in filenames:
    landsat_tif = Path('landsat_scenes/l8_vancouver') / Path(the_file)
    download(str(landsat_tif),dest_folder=dest_folder)
    
band4_bigfile=list(dest_folder.glob("*20180131*B4.TIF"))[0]
band3_bigfile=list(dest_folder.glob("*20180131*B3.TIF"))[0]
meta_bigfile=list(dest_folder.glob("*20180131*MTL.txt"))[0]
with rasterio.open(band4_bigfile) as raster:
    big_transform=raster.transform
    big_profile=raster.profile<jupyter_output>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\Pearl\Desktop\ATSC 301\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_B4.TIF

LC08_L1TP_047026_20150614_20180131_01_T1_B4.TIF already exists
and is 75252265 bytes
will not overwrite

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\Pearl\Desktop\ATSC 301\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131_01_T1_B3.TIF

LC08_L1TP_047026_20150614_20180131_01_T1_B3.TIF already exists
and is 74130783 bytes
will not overwrite

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\Pearl\Desktop\ATSC 301\a301_code\data\landsat8\vancouver\LC08_L1TP_047026_20150614_20180131[...]<jupyter_text>Rename the band 5 tiff file created by the image_zoom notebookChange it from small_file.tiff to the more descriptive band5_reflect.tiffNote that this will fail the second time this cell is executed, so ignore that error<jupyter_code>try:
    tiff_filename=a301.data_dir / Path('small_file.tiff')   
    band5_smallfile=a301.data_dir / Path('band5_reflect.tiff') 
    os.rename(tiff_filename,band5_smallfile)
except FileNotFoundError:
    pass<jupyter_output><empty_output><jupyter_text>Open the small file and get the transform<jupyter_code>with rasterio.open(band5_smallfile) as raster:
    small_transform=raster.transform
    crs=raster.crs
    small_profile=raster.profile
    b5_refl=raster.read(1)<jupyter_output><empty_output><jupyter_text>Use the transform to make a Window with the correct ul corner<jupyter_code>metadata=landsat_metadata(meta_bigfile)    
small_ul_xy = small_transform*(0,0)
big_ul_col_row = ~big_transform*(small_ul_xy[:])
big_col, big_row = int(big_ul_col_row[0]),int(big_ul_col_row[1])
small_height,small_width = b5_refl.shape
small_window=Window(big_col, big_row, small_width, small_height)<jupyter_output><empty_output><jupyter_text>Read only the window pixels from the band 3, 4 files<jupyter_code>refl_dict=dict()
refl_dict[5]=b5_refl
for bandnum,filepath in zip([3,4],[band3_bigfile,band4_bigfile]):
    with rasterio.open(filepath) as src:
        counts = src.read(1, window=small_window)
        refl = calc_reflc_8(counts,bandnum,metadata)
        refl_dict[bandnum]=refl<jupyter_output><empty_output><jupyter_text>check to make sure the read worked<jupyter_code>vmin=0.
vmax=0.1
cmap_ref=plt.get_cmap('viridis')
cmap_ref.set_over('c')
cmap_ref.set_under('b',alpha=0.2)
cmap_ref.set_bad('0.75') #75% grey
the_norm=Normalize(vmin=vmin,vmax=vmax,clip=False)
fig, ax1 = plt.subplots(1,1,figsize=(15,10))
cax=ax1.imshow(refl_dict[3],norm=the_norm,cmap=cmap_ref)
fig.colorbar(cax,ax=ax1)
ax1.set_title('band 3 green reflectance');<jupyter_output><empty_output><jupyter_text>Compare the green and red bands<jupyter_code>fig, (ax1,ax2,ax3) = plt.subplots(1,3,figsize=(15,4))
diff = refl_dict[5] - refl_dict[4]
ax1.hist(refl_dict[3][~np.isnan(refl_dict[3])].flat)
ax2.hist(refl_dict[4][~np.isnan(refl_dict[4])].flat)
ax3.hist(diff[~np.isnan(diff)].flat)
ax1.set(title='band 3 green refl',xlim=[0,0.3])
ax2.set(title='band 4 red refl',xlim=[0,0.3])
ax3.set(title='red - green',xlim=[-0.2,0.2]);<jupyter_output><empty_output><jupyter_text>Write out a new compositeESRI suggests that vegetation looks good if RGB is mapped to [bands 5,4,3](https://www.esri.com/arcgis-blog/products/product/imagery/band-combinations-for-landsat-8).  Reuse the code from the rasterio_demo notebook<jupyter_code>from skimage import  img_as_ubyte
from skimage import exposure

height,width = refl_dict[5].shape
num_chans=3
channels=np.empty([num_chans,height,width],dtype=refl_dict[5].dtype)
for index,the_chan in enumerate([5,4,3]):
    channels[index,:,:] = refl_dict[the_chan]

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

        
with rasterio.open(tif_filename) as infile:
    print(f"\nnew profile: {pprint.pformat(infile.profile)}\n")
    profile=infile.profile

int_channels=np.empty([num_chans,height,width],dtype=np.uint8)    
profile['driver']='PNG'  
profile['dtype'] = 'uint8'
png_filename=tif_filename.with_suffix('.png')
with rasterio.open(png_filename, 'w', **profile) as dst:
    for index,chan in enumerate([5,4,3]):
        image = refl_dict[chan]
        stretched=exposure.equalize_hist(image)   
        int_channels[index,:,:]=img_as_ubyte(stretched)[:,:]
    dst.write(int_channels)
Image(str(png_filename))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11162600517922856, 'fraction_numerical': 0.0734632683658171, 'mean_word_length': 6.683769633507853, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '244165', 'license': 'MIT', 'n_tokens_mistral': 3107, 'n_tokens_neox': 2806, 'n_words': 482, 'path': 'notebooks/hires_multichan-vancouver.ipynb', 'repo_name': 'Pearl-Ayem/ATSC_Course_Work', 'size': '7337'}