Прошу прощения, если название неясно, я новичок в python, и мой словарный запас ограничен.
Я пытаюсь применить растяжку стандартного отклонения к каждой группе в .tif растр, а затем создать новый растр (.tif), накладывая эти полосы с помощью GDAL (Python).
Я могу создавать новые растры ложных цветов с различными комбинациями полос и сохранять их, и я могу создайте желаемое ИЗОБРАЖЕНИЕ в python, используя dstack (первый блок кода), но я не могу сохранить это изображение в виде файла с географической привязкой .tif.
Поэтому для создания растянутого изображения с использованием dstack мой код выглядит следующим образом :
import os
import numpy as np
import matplotlib.pyplot as plt
import math
from osgeo import gdal
# code from my prof
def std_stretch_data(data, n=2):
"""Applies an n-standard deviation stretch to data."""
# Get the mean and n standard deviations.
mean, d = data.mean(), data.std() * n
# Calculate new min and max as integers. Make sure the min isn't
# smaller than the real min value, and the max isn't larger than
# the real max value.
new_min = math.floor(max(mean - d, data.min()))
new_max = math.ceil(min(mean + d, data.max()))
# Convert any values smaller than new_min to new_min, and any
# values larger than new_max to new_max.
data = np.clip(data, new_min, new_max)
# Scale the data.
data = (data - data.min()) / (new_max - new_min)
return data
# open the raster
img = gdal.Open(r'/Users/Rebekah/ThesisData/TestImages/OG/OG_1234.tif')
#open the bands
red = img.GetRasterBand(1).ReadAsArray()
green = img.GetRasterBand(2).ReadAsArray()
blue = img.GetRasterBand(3).ReadAsArray()
# create alpha band where a 0 indicates a transparent pixel and 1 is a opaque pixel
# (this is from class and i dont FULLY understand it)
alpha = np.where(red + green + blue ==0, 0, 1).astype(np.byte)
red_stretched = std_stretch_data(red, 1)
green_stretched = std_stretch_data(green, 1)
blue_stretched = std_stretch_data(blue, 1)
data_stretched = np.dstack((red_stretched, green_stretched, blue_stretched, alpha))
plt.imshow(data_stretched)
plt.show()
И это дает мне прекрасное изображение именно того, что я хочу в отдельном окне. Но нигде в этом коде нет возможности назначать проекции или сохранять их как многополосный TIF.
Поэтому я взял это и применил его как можно лучше к коду, который я использую для создания изображений в ложных цветах, и это не удается (код ниже). Если я создаю 4-полосный TIF с альфа-диапазоном, на выходе получается пустой TIF, а если я создаю 3-полосный TIF и опускаю альфа-диапазон, на выходе получается полностью черный TIF.
import os
import numpy as np
import matplotlib.pyplot as plt
import math
from osgeo import gdal
#code from my professor
def std_stretch_data(data, n=2):
"""Applies an n-standard deviation stretch to data."""
# Get the mean and n standard deviations.
mean, d = data.mean(), data.std() * n
# Calculate new min and max as integers. Make sure the min isn't
# smaller than the real min value, and the max isn't larger than
# the real max value.
new_min = math.floor(max(mean - d, data.min()))
new_max = math.ceil(min(mean + d, data.max()))
# Convert any values smaller than new_min to new_min, and any
# values larger than new_max to new_max.
data = np.clip(data, new_min, new_max)
# Scale the data.
data = (data - data.min()) / (new_max - new_min)
return data
#open image
img = gdal.Open(r'/Users/Rebekah/ThesisData/TestImages/OG/OG_1234.tif')
# get geotill driver
gtiff_driver = gdal.GetDriverByName('GTiff')
# read in bands
red = img.GetRasterBand(1).ReadAsArray()
green = img.GetRasterBand(2).ReadAsArray()
blue = img.GetRasterBand(3).ReadAsArray()
# create alpha band where a 0 indicates a transparent pixel and 1 is a opaque pixel
# (this is from class and i dont FULLY understand it)
alpha = np.where(red + green + blue ==0, 0, 1).astype(np.byte)
# apply the 1 standard deviation stretch
red_stretched = std_stretch_data(red, 1)
green_stretched = std_stretch_data(green, 1)
blue_stretched = std_stretch_data(blue, 1)
# create empty tif file
NewImg = gtiff_driver.Create('/Users/riemann/ThesisData/TestImages/FCI_tests/1234_devst1.tif', img.RasterXSize, img.RasterYSize, 4, gdal.GDT_Byte)
if NewImg is None:
raise IOerror('could not create new raster')
# set the projection and geo transform of the new raster to be the same as the original
NewImg.SetProjection(img.GetProjection())
NewImg.SetGeoTransform(img.GetGeoTransform())
# write new bands to the new raster
band1 = NewImg.GetRasterBand(1)
band1.WriteArray(red_stretched)
band2 = NewImg.GetRasterBand(2)
band2.WriteArray(green_stretched)
band3= NewImg.GetRasterBand(3)
band3.WriteArray(blue_stretched)
alpha_band = NewImg.GetRasterBand(4)
alpha_band.WriteArray(alpha)
del band1, band2, band3, img, alpha_band
Я не совсем уверен, как go отсюда и создать новый файл, отображающий растяжение на разных полосах.
Изображение представляет собой просто 4-полосный растр (NAIP), загруженный с earthexplorer, я могу при необходимости загрузить указанное c изображение, которое я использую для своего теста, но в этом файле нет ничего особенного по сравнению с другими NAIP изображения.