Landsat Collection 2: Stacking and subsetting 

Steps:

1. Download data as .tar.gz files

2. Extract so each data is a folder

3. In a UNIX shell with gdal installed, navigate to the folder containing all the data folder and execute the following command:

• Landsat 4-7

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

• Landsat 8-9

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

• Landsat 8-9 C2L2 Bulk Download

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]

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

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

5. In Python, run the following script to convert to [0,1] reflectance values: ds_scaleC2L2.py

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

• Landsat 4-7

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

• Landsat 8-9

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):

1. Make a new conda environment with the following packages installed:

   • gdal, pandas, rasterio

2. Download data as .tar.gz files

3. Extract so each data is a folder

4. In a Windows PowerShell with an anaconda environment activated with gdal installed, navigate to the folder containing all the data folder and execute the following command:

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") }

• Replace the text in red with the path to your own Anaconda environment

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"}

• Replace the text in red with the upper left and lower right corners of your study area

• Option if you want to use lat/lon (in WGS-84 decimal degrees) instead of UTM:

foreach ($i in ls */*_e){gdal_translate.exe -of ENVI -projwin  18.936430 -33.534807 19.343072 -34.142850 -projwin_srs epsg:4326  ${i} "${i}_sub"}

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

• cd to location where you want new environment

• conda create --prefix rasterio_env

• conda activate C:\Users\dan.sousa\Desktop\rasterio_env

• conda install -c conda-forge mamba

• mamba install -c conda-forge rasterio

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}

(I need to fix stuff in red still)

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