<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Introduction2&nbsp;&nbsp;Get bands 3, 4, 5 fullsize (green, red, near-ir)3&nbsp;&nbsp;This cell reads in your affine transform, metadata and profile4&nbsp;&nbsp;This cell gets the right reflection function for your satellite5&nbsp;&nbsp;Read only the window pixels from the band 3, 4 files6&nbsp;&nbsp;In the next cell plot a mapped ndvi image with a red dot in your ul corner and a white dot in your lr corner  IntroductionThere are 4 cells that ask for changes below, the rest should run as long as youuse the variable names I ask for in the questions.<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
from a301.landsat.toa_reflectance import calc_refl_457, calc_reflc_8
from mpl_toolkits.axes_grid1 import make_axes_locatable

#good cell

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 import windows
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
from a301.landsat.toa_reflectance import calc_refl_457, calc_reflc_8
from mpl_toolkits.axes_grid1 import make_axes_locatable<jupyter_output><empty_output><jupyter_text>Get bands 3, 4, 5 fullsize (green, red, near-ir)At the end of this cell you shiould have the following path objects for your spring scene:meta_bigfile, band3_bigfile, band4_bigfile, band5_bigfilethat point to your landsat TIF and mtl.txt files.<jupyter_code>#good cell

filenames=["LC08_L1TP_190031_20170528_20170615_01_T1_B3.TIF",
    "LC08_L1TP_190031_20170528_20170615_01_T1_B4.TIF",
    "LC08_L1TP_190031_20170528_20170615_01_T1_B5.TIF",
    "LC08_L1TP_190031_20170528_20170615_01_T1_MTL.txt"]
dest_folder=a301.data_dir / Path("landsat8/italy")

band3_bigfile=list(dest_folder.glob("*_B3.TIF"))[0]
band4_bigfile=list(dest_folder.glob("*_B4.TIF"))[0]
band5_bigfile=list(dest_folder.glob("*_B5.TIF"))[0]
meta_bigfile=list(dest_folder.glob("*MTL.txt"))[0]<jupyter_output><empty_output><jupyter_text>This cell reads in your affine transform, metadata and profileUsing band4_bigfile (arbitrary)<jupyter_code>metadata=landsat_metadata(meta_bigfile)
with rasterio.open(band4_bigfile) as raster:
    big_transform=raster.transform
    big_profile=raster.profile

zone = metadata.UTM_ZONE  
crs = cartopy.crs.UTM(zone, southern_hemisphere=False)
p_utm=Proj(crs.proj4_init)
p_lonlat=Proj(proj='latlong',datum='WGS84')
#good cell
metadata=landsat_metadata(meta_bigfile)
with rasterio.open(str(band4_bigfile)) as raster:
    big_transform=raster.affine
    big_profile=raster.profile

zone = metadata.UTM_ZONE  
crs = cartopy.crs.UTM(zone, southern_hemisphere=False)
p_utm=Proj(crs.proj4_init)
p_lonlat=Proj(proj='latlong',datum='WGS84')<jupyter_output><empty_output><jupyter_text>This cell gets the right reflection function for your satellite<jupyter_code>refl_dict={'LANDSAT_7':calc_refl_457,'LANDSAT_8':calc_reflc_8}         
satellite=metadata.SPACECRAFT_ID
refl_fun=refl_dict[satellite]<jupyter_output><empty_output><jupyter_text>Define a subscene window and a transformIn the cell below, get the upper left col,row (ul_col,ul_row) and upper left and lowerright x,y (ul_x,ul_y,lr_x,lr_y)coordinates the upper left corner of your subscene as in the image_zoom notebook.  Use ul_col, ul_row, ul_x, ul_y plus your subscenewidth and height to make a rasterio window and new transform.    window=Window(ul_col, ul_row, small_width, small_height)    new_affine=Affine(30.,0.,ul_x,0.,-30.,ul_y)    extent = [ul_x,lr_x,lr_y,ul_y]<jupyter_code>italy_lon = 13.66477
italy_lat = 41.75983
italy_x, italy_y =proj_transform(p_lonlat,p_utm,italy_lon, italy_lat) 


full_ul_xy=np.array(big_transform*(0,0))
print(f"orig ul corner x,y (km)={full_ul_xy*1.e-3}")

ul_col, ul_row = ~big_transform*(italy_x,italy_y)
ul_col, ul_row = int(ul_col), int(ul_row)

l_col_offset= -1300
r_col_offset= +2000
b_row_offset= +2600
t_row_offset= -100
col_slice=slice(ul_col+l_col_offset,ul_col+r_col_offset)
row_slice=slice(ul_row + t_row_offset, ul_row + b_row_offset)
italy_ul_xy = big_transform*(col_slice.start,row_slice.start)
italy_lr_xy = big_transform*(col_slice.stop,row_slice.stop)

small_height, small_width = 2700,3300
ul_x, ul_y = italy_ul_xy[0], italy_ul_xy[1]

# window=Window(ul_col, ul_row, small_width, small_height)
new_affine=Affine(30.,0.,ul_x,0.,-30.,ul_y)
image_extent=[italy_ul_xy[0],italy_lr_xy[0],italy_ul_xy[1],italy_lr_xy[1]]<jupyter_output>orig ul corner x,y (km)=[ 271.185 4742.715]<jupyter_text>Read only the window pixels from the band 3, 4, 5 files<jupyter_code>refl_dict=dict()
for bandnum,filepath in zip([3,4,5],[band3_bigfile,band4_bigfile,band5_bigfile]):
    with rasterio.open(filepath) as src:
        counts = src.read(1, window=small_window)
        refl_vals = refl_fun(counts,bandnum,metadata)
        refl_dict[bandnum]=refl_vals
#good cell
from a301.landsat.toa_reflectance import toa_reflectance_8
refl_vals=toa_reflectance_8([3,4,5],meta_bigfile)
refl_dict=dict()
for bandnum,filepath in zip([3,4,5],[band3_bigfile,band4_bigfile,band5_bigfile]):
    with rasterio.open(str(filepath)) as src:
#         counts = src.read(1, window=small_window)
#         refl_vals = (,bandnum,metadata)
        refl_dict[bandnum]=refl_vals<jupyter_output>Scene LC81900312017148LGN00 center time is 2017-05-28 09:46:46<jupyter_text>In the next cell calculate your ndviSave it in a variable called ndvi<jupyter_code># YOUR CODE HERE
ndvi = (refl_vals[5] - refl_vals[4])/(refl_vals[5] + refl_vals[4])
plt.hist(ndvi[~np.isnan(ndvi)].flat);
plt.title('spring ndvi')
plt.savefig('spring_ndvi.png')<jupyter_output><empty_output><jupyter_text>In the next cell plot a mapped ndvi image with a red dot in your ul corner and a white dot in your lr cornerAdjust this plot to fit your image.  Just delete the bottom line and work with the provided commands<jupyter_code>vmin=0.0
vmax=0.8
the_norm=Normalize(vmin=vmin,vmax=vmax,clip=False)
palette='viridis'
pal = plt.get_cmap(palette)
pal.set_bad('0.75') #75% grey for out-of-map cells
pal.set_over('w')  #color cells > vmax red
pal.set_under('k')  #color cells < vmin black
fig, ax = plt.subplots(1, 1,figsize=[10,15],
                       subplot_kw={'projection': crs})
col=ax.imshow(ndvi,origin="upper",
         extent=image_extent,transform=crs)
ax.plot(italy_ul_xy[0],italy_ul_xy[1],'ro',markersize=50)
ax.plot(italy_lr_xy[0],italy_lr_xy[1],'wo',markersize=50)
ax.set(title="spring ndvi")
cbar_ax = fig.add_axes([0.95, 0.2, 0.05, 0.6])
cbar=ax.figure.colorbar(col,extend='both',cax=cbar_ax,orientation='vertical')
cbar.set_label('ndvi index')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1008174386920981, 'fraction_numerical': 0.05013623978201635, 'mean_word_length': 7.093715545755237, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '82066', 'license': 'MIT', 'n_tokens_mistral': 3011, 'n_tokens_neox': 2792, 'n_words': 598, 'path': 'assignments/ndvi_subscene_assignment-PA_v1.ipynb', 'repo_name': 'Pearl-Ayem/ATSC_Course_Work', 'size': '7340'}