I am trying to calculate the dNBR of a particular area using Landsat 8 in python and save the output as a GeoTiff so I can use it in QGIS.
Following other tutorials the dNBR calculation works fine but I am struggling to now save this variable (dnbr_landsat) as a GeoTiff. It is creating a GeoTiff but when I view it in QGIS it's completely black with no values like it hasn't actually been assigned the array, and in the console I get the following error on Line 103. "AttributeError: 'Band' object has no attribute 'writeArray'"
os.chdir("C:/Users/my_name/Landsat_08")
# Stack the Landsat 8 bands
# This creates a numpy array with each "layer" representing a single band
landsat_prefire_path = glob(
"LC08_L1TP_095086_20141221_20200910_02_T1/LC08_L1TP_095086_20141221_20200910_02_T1*.tif"
)
prefire_bands = landsat_prefire_path
# We handle the connections with "with"
with rasterio.open(prefire_bands[4]) as src:
prefire_b5 = src.read(1)
with rasterio.open(prefire_bands[6]) as src:
prefire_b7 = src.read(1)
# Allow division by zero
np.seterr(divide='ignore', invalid='ignore')
# Calculate NDVI
prefire_nbr = (prefire_b5.astype(float) - prefire_b7.astype(float)) / (prefire_b5 + prefire_b7)
# Take a spatial subset of the ndvi layer produced
ndvi_sub = prefire_nbr[2000:3000, 2000:3000]
# Plot
plt.imshow(ndvi_sub)
plt.show()
landsat_postfire_path = glob("LC08_L1TP_095086_20150106_20200910_02_T1/LC08_L1TP_095086_20150106_20200910_02_T1*.tif")
data = gdal.Open("C:/Users/my_name/Landsat_08/LC08_L1TP_095086_20150106_20200910_02_T1/LC08_L1TP_095086_20150106_20200910_02_T1_b5.tif")
band = data.GetRasterBand(1)
arr = band.ReadAsArray()
[cols, rows] = arr.shape
postfire_bands = landsat_postfire_path
with rasterio.open(postfire_bands[4]) as src:
postfire_b5 = src.read(1)
with rasterio.open(postfire_bands[6]) as src:
postfire_b7 = src.read(1)
# Allow division by zero
np.seterr(divide='ignore', invalid='ignore')
# Calculate NDVI
postfire_nbr = (postfire_b5.astype(float) - postfire_b7.astype(float)) / (postfire_b5 + postfire_b7)
# Take a spatial subset of the ndvi layer produced
ndvi_sub = postfire_nbr[2000:3000, 2000:3000]
# Plot
plt.imshow(ndvi_sub)
plt.show()
dnbr_landsat = prefire_nbr - postfire_nbr
# Define dNBR classification bins
dnbr_class_bins = [-np.inf, -.1, .1, .27, .66, np.inf]
#dnbr_landsat_class = np.digitize(dnbr_landsat, dnbr_class_bins)
dnbr_landsat_class = xr.apply_ufunc(np.digitize,
dnbr_landsat,
dnbr_class_bins)
plt.imshow(dnbr_landsat_class)
plt.show()
gt = data.GetGeoTransform()
proj = data.GetProjection()
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create("dNBR_Landsat2.tif", rows, cols, 1, gdal.GDT_UInt16)
dataset.SetGeoTransform(gt)
dataset.SetProjection(proj)
tif_data = dataset.GetRasterBand(1)
tif_data.writeArray(dnbr_landsat)
tif_data.SetNoDataValue(np.nan)
tif_data.FlushCache()
tif_data = None
dataset = None
Relatively new to Python, so I'm sure it might be something simple I've missed.