Landsat Collection 2: Stacking and subsetting 

Steps:

for i in */*LT0[4-5]*B1.TIF; do gdal_merge.py -of ENVI -o ${i%???????}_e ${i} ${i%?????}2.TIF ${i%?????}3.TIF ${i%?????}4.TIF ${i%?????}5.TIF ${i%?????}7.TIF  -separate; done

for i in */*LC0[8-9]*B1.TIF; do gdal_merge.py -of ENVI -o ${i%???????}_e ${i} ${i%?????}2.TIF ${i%?????}3.TIF ${i%?????}4.TIF ${i%?????}5.TIF ${i%?????}6.TIF ${i%?????}7.TIF  -separate; done

VSWIR - [for i in *LC0[8-9]*B2.TIF; do gdal_merge.py -of ENVI -o ${i%???????}_e ${i} ${i%?????}3.TIF ${i%?????}4.TIF ${i%?????}5.TIF ${i%?????}6.TIF ${i%?????}7.TIF  -separate; done

Thermal - [for i in *LC0[8-9]*B10.TIF; do gdal_merge.py -of ENVI -o ${i%???????}e_T ${i} -separate; done]


for i in */*_e; do gdal_translate -of ENVI -projwin 210000 3970000 350000 3880000 ${i} ${i}_sub; done

for i in */*_e.hdr; do sed 's/data type = 12/data type = 4/' ${i} > ${i%????}Scale.hdr; done

for i in */*Scale.hdr; do echo 'wavelength units = Micrometers' >> ${i}; done

for i in */*Scale.hdr; do echo 'wavelength = {0.48, 0.56, 0.665, 0.83, 1.6, 2.2}' >> ${i}; done

for i in */*Scale.hdr; do echo 'wavelength = {0.44, 0.48, 0.56, 0.665, 0.83, 1.6, 2.2}' >> ${i}; done

Band centers from: https://www.usgs.gov/faqs/what-are-band-designations-landsat-satellites 

New instructions for Windows machines in Storm Hall (Feb 7 2024):


foreach ($i in ls */*B1.TIF){python C:\Users\dan.sousa\Desktop\gdal_env\Scripts\gdal_merge.py -o  ($i.DirectoryName + '\' + $i.Name).Replace("\","\\").Replace("B1.TIF","_e") -of ENVI -separate ($i.DirectoryName + '\' + $i.Name).Replace("\","\\").Replace("B1.TIF","B1.TIF") ($i.DirectoryName + '\' + $i.Name).Replace("\","\\").Replace("B1.TIF","B2.TIF") ($i.DirectoryName + '\' + $i.Name).Replace("\","\\").Replace("B1.TIF","B3.TIF") ($i.DirectoryName + '\' + $i.Name).Replace("\","\\").Replace("B1.TIF","B4.TIF") ($i.DirectoryName + '\' + $i.Name).Replace("\","\\").Replace("B1.TIF","B5.TIF") ($i.DirectoryName + '\' + $i.Name).Replace("\","\\").Replace("B1.TIF","B6.TIF") ($i.DirectoryName + '\' + $i.Name).Replace("\","\\").Replace("B1.TIF","B7.TIF") }


4.(Optional) At the same location execute the following command, replacing the numbers with bounding X and Y coordinates of your spatial subset: 

foreach ($i in ls */*_e){gdal_translate.exe -of ENVI -projwin 470000 3636000 490000 3623000 ${i} "${i}_sub"}


5. Make a new environment with 'rasterio' installed; activate it


6. Download the following script to convert to [0,1] reflectance values: ds_scaleC2L2_lab.py. Run it:

Change the line that says directory = 'C:\\Users\\dan.sousa\\Downloads\\'  to the directory that has your own data.


7. Now run the following

foreach ($i in ls */*_e_sub.hdr){cp -v $i ($i.DirectoryName + '\' + $i.Name).Replace('sub.hdr','subScale.hdr')}

foreach ($i in ls */*_e_subScale.hdr){(Get-Content -Raw $i) -replace 'data type = 12', 'data type = 4' | Set-Content -NoNewLine $i}


for i in */*Scale.hdr; do echo 'wavelength units = Micrometers' >> ${i}; done

for i in */*Scale.hdr; do echo 'wavelength = {0.44, 0.48, 0.56, 0.665, 0.83, 1.6, 2.2}' >> ${i}; done