0

I'm trying to determine if my code for determining land surface temperature is correct. The imagery contains a fire in the middle of October with snow mixed in so temperatures show be nice and contrasting but when I go to view the created image in QGIS, the image is showing up as black. I realize that this doesn't mean there isn't anything there, I just want to make sure it's the code that is or is not the issue.

import os
import rasterio
import glob
from rasterio.merge import merge
import matplotlib.pyplot as plt
from rasterio.plot import show
import numpy as np

# Directory
# These directories are listed as "A" and "B" because one image doesn't cover the 
# area of interest. So there are two sets and we will mosaic them later. "A" and "B" 
# don't mean anything more than to organize them and show they are two different datasets



# Band Specifications - Surface Reflectance
# Band 1 = Coastal Aerosol 
# Band 2 = Blue
# Band 3 = Green
# Band 4 = Red 
# Band 5 = NIR
# Band 6 = SWIR 1 (1.57 - 1.65)
# Band 7 = SWIR 2 (2.11 - 2.29)

# Thermal Band
# Band 10 = TIRS 1 (10.6 - 11.19)


# Lists and Bands
# I decided to read in all of the bands from the datasets before the scripts that create 
# our products for analysis before the section that would create those products.  This 
# was done so that access to all the data is in one place rather than sporadically located
# throughout the script.  This way also ensures common and organized variable naming 
# conventions. 

# Like the "A" and "B" directories, the lists and bands are named to match the directory.

duringA = [
    'LC08_L2SP_035031_20201013_20201105_02_T1_SR_B1.TIF',
    'LC08_L2SP_035031_20201013_20201105_02_T1_SR_B2.TIF',
    'LC08_L2SP_035031_20201013_20201105_02_T1_SR_B3.TIF',
    'LC08_L2SP_035031_20201013_20201105_02_T1_SR_B4.TIF',
    'LC08_L2SP_035031_20201013_20201105_02_T1_SR_B5.TIF',
    'LC08_L2SP_035031_20201013_20201105_02_T1_SR_B6.TIF',
    'LC08_L2SP_035031_20201013_20201105_02_T1_SR_B7.TIF',
    'LC08_L2SP_035031_20201013_20201105_02_T1_ST_B10.TIF'
    ]

duringAband1= rasterio.open(os.path.join(duringDirA, duringA[0]), mode = 'r')
duringAband2= rasterio.open(os.path.join(duringDirA, duringA[1]), mode = 'r') 
duringAband3= rasterio.open(os.path.join(duringDirA, duringA[2]), mode = 'r')
duringAband4= rasterio.open(os.path.join(duringDirA, duringA[3]), mode = 'r')
duringAband5= rasterio.open(os.path.join(duringDirA, duringA[4]), mode = 'r')
duringAband6= rasterio.open(os.path.join(duringDirA, duringA[5]), mode = 'r')
duringAband7= rasterio.open(os.path.join(duringDirA, duringA[6]), mode = 'r')
duringAband10= rasterio.open(os.path.join(duringDirA, duringA[7]), mode = 'r')  
# 1.) TOA

# Radiance Multiplicative rescaling factor.
rmb= 3.3420E-04

#radiance additive rescaling factor
rab = 0.10000

#band 10 'A' top of atmostphere
b10ATOA = (duringAband10.read()*rmb)+rab

# 2.) Conversion of radiance at sensor (BT))
#BT = (k2)/(np.log((k1/TOA)+1))

k1 = 774.8853
k2 = 1321.0789

bT1 = ((k2)/(np.log((k1/b10ATOA)+1))-273.15)

# 3.) Calculate NDVI

# ndvi = ((nir(band5) - red(band4))/(nir(band5) + red(band4)))

#ndvi 
ndvi1 = ((duringAband5.read() - duringAband4.read())/(duringAband5.read() + duringAband4.read()))

# 4.) Calculate the proportion of vegetation (pv))

#pv = square((ndvi - ndvi.min)/(ndvi.max - ndvi.min))
pv1 = np.sqrt(((ndvi1 - np.min(ndvi1))/(np.max(ndvi1) - np.min(ndvi1))))

# 5.) Calculate Land Surface Emissivity
# referred to as 'e' for emissivity
# Equation
# 0.004 * pv + 0.986

# emissivity 1
e1 = 0.004 * pv1 + 0.986

# Land Surface Temperature
lst1 = (bT1 / (1 + (0.00115 * bT1 / 1.4388) * np.log(e1)))

#creating the tifs
#LST 1 stack
lst1file = rasterio.open(
    os.path.join(duringTempDir, 'lst1file.tif'), mode = 'w',
    driver = 'GTiff',
    height = '7981',
    width = '7871',
    count = 1,
    dtype = 'float64',
    crs = duringAband10.crs,
    transform=duringAband10.transform   
    )

lst1file.write(lst1[0], 1)
lst1file.close

lst1file = rasterio.open(os.path.join(duringTempDir, 'lst1file.tif'), mode = 'r')
1
  • For code review there is the Code Review Stack Exchange but here you'll need to we ask that you show us where you are stuck with your own code by including just a code snippet in your question. At the moment we cannot see where duringAband10 is defined.
    – PolyGeo
    Aug 29, 2022 at 1:20

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.