Python 3：如何使用颜色表模糊GeoTIFF图像？

Question

我有一个GeoTIFF图像，我需要通过应用平滑滤镜来模糊。图像本身包含需要保留的元数据。它的位深度为8，并使用具有256个32位RGBA值的颜色表来查找每个像素的颜色，但为了使得到的图像看起来平滑，它可能必须使用比特深度的24或32并且没有颜色表，或者使用jpeg压缩。可能进一步复杂化的是图像大23,899x18,330像素，这几乎是PIL默认打开的最大文件的五倍。

如何在Python 3中创建此图像的模糊版本？

我也尝试使用PIL打开并再次保存：

from PIL import Image
Image.MAX_IMAGE_PIXELS = 1000000000
im = Image.open(file_in)
im.save(file_out)

这段代码没有崩溃，我得到一个新的.tif文件，大致与原始文件一样大，但是当我尝试在Windows Photo Viewer中打开它来查看它时，应用程序说它已损坏了，它不能由PIL重新开放。

我也尝试过使用GDAL。当我尝试this code时，我得到一个大小为835 MB的输出图像，它对应于位深度为16的未压缩图像（这也是当我右键单击它并选择“属性”时文件元数据所说的内容。 - 我正在使用Windows 10）。但是，生成的图像是单色的，非常暗，颜色看起来像是乱七八糟，这使我相信我正在尝试的代码将像素值解释为强度值而不是表键。

因此，为了使这个方法有效，我需要弄清楚如何将颜色表（对于元组的某种容器，类型为osgeo.gdal.ColorTable）应用于栅格波段（无论栅格波段是什么），这是一个numpy形状为(18330, 23899)的数组，以获得一个形状为(18330, 23899, 4)或(4, 18330, 23899)的新numpy数组（不知道哪个是正确的形状），将其插回到加载的图像中并删除颜色表（或者创建一个新的颜色表相同的元数据），最后保存已启用压缩的修改后的图像（所以我更接近原始文件大小 - 11.9 MB - 而不是835 MB，这是我现在得到的文件的大小）。我怎样才能做到这一点？

Answer 1

pyvips可以使用少量内存快速处理大量图像，并支持调色板TIFF图像。

不幸的是，它不支持额外的geotiff标签，因为libtiff不适用于未知标签类型。您需要以其他方式复制该元数据。

无论如何，如果你能做到这一点，pyvips应该适用于你的形象。我试过这个例子：

import sys
import pyvips

# the 'sequential' hint tells libvips that we want to stream the image
# and don't need full random access to pixels ... in this mode, 
# libvips can read, process and write in parallel, and without needing
# to hold the whole image in memory
image = pyvips.Image.new_from_file(sys.argv[1], access='sequential')
image = image.gaussblur(2)
image.write_to_file(sys.argv[2])

在您拥有的类型和大小的图像上，生成JPEG压缩的TIFF：

$ tiffinfo x2.tif 
TIFF Directory at offset 0x1a1c65c6 (438068678)
  Image Width: 23899 Image Length: 18330
  Resolution: 45118.5, 45118.5 pixels/cm
  Bits/Sample: 8
  Compression Scheme: None
  Photometric Interpretation: palette color (RGB from colormap)
...
$ /usr/bin/time -f %M:%e python3  ~/try/blur.py x2.tif x3.tif[compression=jpeg]
137500:2.42

所以140MB的内存，2.5秒。输出图像看起来是正确的，是24mb，所以离你的距离太近了。

Answer 2

栅格波段只是图像的每个“图层”的名称，在您的情况下，它们将是红色，绿色，蓝色和Alpha值。这些是你想要模糊的东西。您可以使用data.GetRasterBand(i)打开图像并将每个波段保存到单独的数组中，以获得使用GDAL打开的图像的第i个波段（使用1索引，而不是0索引）。

然后，您可以尝试使用SciPy的scipy.ndimage.gaussian_filter来实现模糊。您需要向它发送一个形状为（x，y）的数组，因此您必须分别为每个栅格波段执行此操作。您应该能够使用GDAL将数据另存为另一个GeoTIFF。

---

如果您正在使用的颜色表意味着您的数据以某种奇怪的格式存储在每个栅格波段中，而不是仅为R，G，B和A中的每一个在0和1之间浮动，那么请考虑使用scipy.ndimage.generic_filter，尽管在不知道您的数据如何存储的情况下，很难详细说明您的数据。