Сохранение карты в TIF-файле

[bgnik]

Уважаемые разработчики! Очень часто при формировании растровой подложки приходится использовать отсканированное изображение топографического плана. Обычно этот файл выполняется в ч/б варианте в формате TIF и сжатии методом Group4 Fax. Было бы очень полезно иметь возможность сохранять фрагменты карты в этом формате (с помощью image.save), задавая число точек, bbox и метод сжатия.
Как я понимаю, на сегодня формат TIF провайдером GDAL не поддерживается.
Есть ли перспектива решения этой проблемы в обозримом будущем?

[Александр]

Ниже приведен пример того, как сохранить текущее окно карты в формате TIFF

Code: Select all

from pathlib import Path
from typing import cast

import axipy as axp
from osgeo import gdal
import numpy as np

import PySide2.QtGui as Qg
from PySide2.QtCore import Qt


def get_mapview_image(mv: axp.MapView, dpi: float) -> Qg.QImage:
k = dpi / 96
img = Qg.QImage(int(mv.device_rect.width * k), int(mv.device_rect.height * k), Qg.QImage.Format_RGB32)
img.fill(Qt.white)
context = axp.Context(Qg.QPainter(img))
context.coordsystem = mv.coordsystem
context.rect = mv.scene_rect
context.dpi = dpi
mv.map.draw(context)

return img


def save_image_to_mono_geotiff(img: Qg.QImage, filename: Path, dpi: float) -> None:
driver: gdal.Driver = gdal.GetDriverByName("GTiff")
height = img.height()
width = img.width()
bands = 1

data = cast(memoryview, img.constBits()).tobytes()
np_data = np.frombuffer(data, dtype=np.uint8)
np_data = np_data.reshape(4, width, height, order='F').transpose(0, 2, 1)
np_data = np.mean([np_data[0], np_data[1], np_data[3]], axis=0)
np_data = np.where(np_data == 255, 1, 0)

try:
ds: gdal.Dataset = driver.Create(
str(filename),
width,
height,
bands,
gdal.GDT_Byte,
options=['COMPRESS=CCITTFAX4', 'NBITS=1'],
)
ds.SetMetadata({'TIFFTAG_XRESOLUTION': dpi, 'TIFFTAG_YRESOLUTION': dpi})
ds.WriteArray(np_data)
finally:
ds = None # noqa


# Пример использования.
# Сохраняем текущее окно карты используя tif Group4 Fax.
image = get_mapview_image(axp.view_manager.active, 300)
save_image_to_mono_geotiff(image, Path(r"c:\tmp\test.tif"), 300)
print("Completed.")


[bgnik]

Большое спасибо за оперативность! Попробую и о результатах напишу.

[bgnik]

Проверил, работает. Понравилось быстродействие. То, что я реализовывал на ImageMagic или с помощью PIL работает медленнее. Но есть общая проблема, которую я не могу преодолеть: в ряде программ созданный TIF открывается в инверсном виде. Это, в частности, Microsoft office Picture Manager (это бы ладно, но вот Microstation Descartes - тут проблема). Где то в файле есть признак, что считать черным цветом, 0 или 1. Как до него добраться я не смог понять. Нет ли такой настройки в параметрах провайдера?

[Александр]

Можно попробовать рисовать инвертировано и установить параметр 'PHOTOMETRIC=MINISWHITE'.
Ниже пример конвертации. Кроме того, я внес некоторые улучшения в расчеты для повышения скорости выполнения, если это важно.
Результат я не тестировал. У меня и предыдущий результат открывается нормально.

Code: Select all

from pathlib import Path
from typing import cast

import axipy as axp
from osgeo import gdal
import numpy as np

import PySide2.QtGui as Qg
from PySide2.QtCore import Qt


def get_mapview_image(mv: axp.MapView, dpi: float) -> Qg.QImage:
k = dpi / 96
img = Qg.QImage(int(mv.device_rect.width * k), int(mv.device_rect.height * k), Qg.QImage.Format_RGB32)
img.fill(Qt.white)
context = axp.Context(Qg.QPainter(img))
context.coordsystem = mv.coordsystem
context.rect = mv.scene_rect
context.dpi = dpi
mv.map.draw(context)

return img


def save_image_to_mono_geotiff(img: Qg.QImage, filename: Path, dpi: float) -> None:
driver: gdal.Driver = gdal.GetDriverByName("GTiff")
height = img.height()
width = img.width()
bands = img.depth() // 8

data = cast(memoryview, img.constBits()).tobytes()
np_data = np.frombuffer(data, dtype=np.uint8)
np_data = 1 - np.min(np_data.reshape(height, width, bands), 2) // 255
try:
ds: gdal.Dataset = driver.Create(
str(filename),
width,
height,
1,
gdal.GDT_Byte,
options=['COMPRESS=CCITTFAX4', 'NBITS=1', 'PHOTOMETRIC=MINISWHITE']
)
ds.SetMetadata({'TIFFTAG_XRESOLUTION': dpi, 'TIFFTAG_YRESOLUTION': dpi})

ds.WriteArray(np_data)
ds.FlushCache()
finally:
ds = None # noqa


# Пример использования.
# Сохраняем текущее окно карты используя tif Group4 Fax.
image = get_mapview_image(axp.view_manager.active, 300)
save_image_to_mono_geotiff(image, Path(r"c:\tmp\test.tif"), 300)
print("Completed.")



[bgnik]

Вот теперь сохраняется нормально! И в аксиоме, и в FastStone, и в Picture Manager, и в Microstation открывается одинаково правильно. Предыдущий вариант открывался в аксиоме и в FastStone нормально, а в Picture Manager и в Microstation - инверсно.
Быстродействие отличное! Теперь у меня другая проблема: мне нужно сохранить не окно, а фрагмент карты, полученный с помощью image=map.to_image(width, height, coordsystem,bbox). Напрямую засунуть эту image в save_image_to_mono_geotiff(image, Path(r"c:\tmp\test.tif"), 300) не получается (наверно, это разные объекты - ?), что делается внутри этой функции понимаю плохо (по-видимому что-то не учитываю с Numpy - слабо её знаю). Подскажите, пожалуйста, и я буду счастлив!

[Александр]

Map.to_image(width, height, coordsystem,bbox) сохраняет картинку на прозрачном фоне, а не на белом.
Функция save_image_to_mono_geotiff рассчитана на то, что белый цвет будет оставаться белым, а остальные станут чёрными.
Соответственно прозрачный цвет также станет чёрным.

Проще всего не использовать Map.to_image(), а использовать Map.draw()

Code: Select all

def get_map_image(map_: axp.Map, width: int, height: int, coordsystem: axp.CoordSystem, bbox: axp.Rect, dpi: float) -> Qg.QImage:
img = Qg.QImage(width, height, Qg.QImage.Format_ARGB32)
img.fill(Qt.white)
context = axp.Context(Qg.QPainter(img))
context.coordsystem = coordsystem
context.rect = bbox
context.dpi = dpi
map_.draw(context)

return img

# Пример использования.
mv = axp.view_manager.active
width = 1024
height = 768
coordsystem = mv.coordsystem
bbox = mv.scene_rect

image = get_map_image(mv.map, width, height, coordsystem, bbox, 96)
save_image_to_mono_geotiff(image, Path(r"c:\tmp\test.tif"), 96)
print("Completed.")



[bgnik]

Отлично! Огромное спасибо! Все как надо: и быстро, и корректно. Одно плохо - что делается в save_image_to_mono_geotiff для меня загадка, использовал втупую. Ничего, разберусь со временем. Впереди еще одна проблема: как мне засунуть в эту функцию не image из axipy, а Image из PIL?

[Александр]

Одно плохо - что делается в save_image_to_mono_geotiff для меня загадка, использовал втупую. Ничего, разберусь со временем.

Прокомментирую эту часть:

Code: Select all

data = cast(memoryview, img.constBits()).tobytes()
np_data = np.frombuffer(data, dtype=np.uint8)
np_data = 1 - np.min(np_data.reshape(height, width, bands), 2) // 255


1. получаем массив байт изображения
- img.constBits() возвращает указатель на данные изображения img.
- cast(memoryview ...) используется для приведения объекта памяти к типу memoryview. Это необязательная операция. Предназначена для того чтобы IDE узнало о возвращаемом типе данных.
- tobytes() преобразует memoryview в байты

про memoryview и bytes можно прочитать здесь:
https://docs.python.org/3/library/stdtypes.html#binary-sequence-types-bytes-bytearray-memoryview

2. создает массив NumPy используя тип данных np.uint8 (беззнаковый 8-битный целочисленный тип).

3. Создаётся новый массив NumPy, где 0 - если исходный пиксел был белым и 1 в противном случае
- np_data.reshape(height, width, bands) изменяет форму массива np_data на размеры (height, width, bands).
- np.min(..., 2) находит минимальное значение по оси 2 (по каналам цвета) массива.
- 1 - ... // 255 выполняет операции для создания нового массива, где пиксели белого цвета становятся 0, а все остальные пиксели становятся 1. У белого цвета минимальным значениям по каналам будет 255. Так что можно поделить все значения массива на 255 с отбрасыванием дробной части и вычесть из единицы. Вместо этого можно сравнить каждый элемент массива с 255 так: != 255, но на моём компьютере это работает медленнее.
Таким образом, каждый пиксель будет иметь значение 0 или 1 в зависимости от того, является ли пиксель белым или нет.

Впереди еще одна проблема: как мне засунуть в эту функцию не image из axipy, а Image из PIL?

Получить NumPy массив из изображения PIL очень просто. Достаточно вызвать numpy.asarray(img)

Code: Select all

import numpy as np
from PIL import Image

img = Image.open(r"c:\tmp\test.tif")
np_data = np.asarray(img)


[bgnik]

Александр, благодарю! Хочу отметить, что ваша команда, в отличие от многих, очень внимательно относится к поднимаемым на форуме вопросам. Это очень ценно, еще раз спасибо за это внимание! Аксиоме это идет только на пользу. Так держать!