Блог

Я пытаюсь рассчитать dNBR определенной области, используя Landsat 8 на python, и сохранить результат в виде геотифа, чтобы я мог использовать его в QGIS.

Следуя другим учебным пособиям, расчет dNBR работает нормально, но я изо всех сил пытаюсь сохранить эту переменную (dnbr_landsat) в качестве геотифа. Он создает геотиф, но когда я просматриваю его в QGIS, он полностью черный, без значений, как будто ему на самом деле не был присвоен массив, и в консоли я получаю следующую ошибку в строке 103. «Ошибка атрибута: объект «Полоса» не имеет атрибута»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  

Относительно новичок в Python, поэтому я уверен, что это может быть что-то простое, что я пропустил.