Задать вопрос
@zedxxx
python

Как работать с растрами в памяти в GDAL (Python 2.x)?

Стоит задача перепроецирования растров на лету (в памяти) средствами GDAL. Т. е. на вход подаётся некое изображение в формате PNG/JPEG/GIF в известной проекции и с известными координатами, а на выходе ожидаем получить изображение в том же формате, но в новой проекции. И вот на моменте передачи растра в GDAL и получении его обратно возникли затруднения. В коде ниже, изображение открывается с диска и в принципе может быть сохранено на диск, но в формате BMP (если dest_ds создавать с соответствующим драйвером), а в PNG или JPEG сохранять не хочет. Если я правильно понимаю, тут нужно использовать какие-то сторонние средства (какие?), чтобы открыть растр, а затем уже использовать gdal.WriteRaster() для передачи его в GDAL. А на выходе, соответственно, нужно считывать растр через метод gdal.ReadRaster() и теми же сторонними средствами паковать его в тот или иной формат? 

import osr
import gdal


def reproject(src_epsg, src_gt, dest_epsg, dest_gt, dest_width, dest_height):

    src_ds = gdal.Open('test.png')

    src_proj = osr.SpatialReference()
    src_proj.ImportFromEPSG(src_epsg)

    src_ds.SetGeoTransform(src_gt)
    src_ds.SetProjection(src_proj.ExportToWkt())

    driver = gdal.GetDriverByName('MEM')
    dest_ds = driver.Create('', dest_width, dest_height, src_ds.RasterCount, gdal.GDT_Byte)

    dest_proj = osr.SpatialReference()
    dest_proj.ImportFromEPSG(dest_epsg)
    
    dest_ds.SetGeoTransform(dest_gt)
    dest_ds.SetProjection(dest_proj.ExportToWkt())

    gdal.ReprojectImage(src_ds, dest_ds, src_proj.ExportToWkt(),
                        dest_proj.ExportToWkt(), gdal.GRA_Bilinear)

    ... ???
  • Вопрос задан
  • 3587 просмотров
Комментировать
Подписаться 2 Оценить Комментировать
Решения вопроса 1
@zedxxx Автор вопроса
Придумал загрузку и выгрузку растров через Pillow и numpy. Возможно, не совсем оптимально с точки зрения производительности, но работает:
warp.py
#!/usr/bin/python
# -*- coding: utf-8 -*-

import osr
import gdal
import numpy
from PIL import Image
from StringIO import StringIO


def reproject_vrt(src_epsg, src_gt, src_size, dst_epsg):

    src_proj = osr.SpatialReference()
    src_proj.ImportFromEPSG(src_epsg)

    driver = gdal.GetDriverByName('VRT')
    src_ds = driver.Create('', src_size[0], src_size[1])

    src_ds.SetGeoTransform(src_gt)
    src_ds.SetProjection(src_proj.ExportToWkt())

    dst_proj = osr.SpatialReference()
    dst_proj.ImportFromEPSG(dst_epsg)

    dst_ds = gdal.AutoCreateWarpedVRT(
        src_ds,
        src_proj.ExportToWkt(),
        dst_proj.ExportToWkt())

    return dst_ds.GetGeoTransform(), dst_ds.RasterXSize, dst_ds.RasterYSize


def reproject_band(band, src_epsg, src_gt, dst_epsg, dst_gt, dst_width, dst_height):

    src_height, src_width = band.shape

    driver = gdal.GetDriverByName('MEM')

    src_ds = driver.Create('', src_width, src_height)
    src_ds.GetRasterBand(1).WriteArray(band)

    src_proj = osr.SpatialReference()
    src_proj.ImportFromEPSG(src_epsg)

    src_ds.SetGeoTransform(src_gt)
    src_ds.SetProjection(src_proj.ExportToWkt())

    dst_ds = driver.Create('', dst_width, dst_height)

    dst_proj = osr.SpatialReference()
    dst_proj.ImportFromEPSG(dst_epsg)

    dst_ds.SetGeoTransform(dst_gt)
    dst_ds.SetProjection(dst_proj.ExportToWkt())

    gdal.ReprojectImage(src_ds, dst_ds, None, None, gdal.GRA_Bilinear)

    return dst_ds.ReadAsArray()


def reproject_raster(raster, src_epsg, src_gt, dst_epsg):

    src_img = Image.open(StringIO(raster))

    dst_gt, dst_width, dst_height = reproject_vrt(src_epsg, src_gt, src_img.size, dst_epsg)

    src_array = numpy.array(src_img)

    if src_array.ndim == 3:
        dst_array_shape = (dst_height, dst_width, src_array.shape[2])
        dst_array = numpy.zeros(dst_array_shape, src_array.dtype)
        for i in range(0, src_array.shape[2]):
            band = src_array[:, :, i]
            band = reproject_band(band, src_epsg, src_gt, dst_epsg, dst_gt, dst_width, dst_height)
            dst_array[:, :, i] = band
    elif src_array.ndim == 2:
        dst_array = reproject_band(src_array, src_epsg, src_gt, dst_epsg, dst_gt, dst_width, dst_height)
    else:
        raise Exception('Unexpected array geometry!')

    dst_img = Image.frombytes(src_img.mode, (dst_width, dst_height), dst_array.data)

    raster = StringIO()
    dst_img.save(raster, format=src_img.format)

    return raster.getvalue()

И ещё, как я понял перепроецировать нужно каждый канал (RGB) отдельно, иначе получается каша...
Ответ написан
Комментировать
Комментировать
Пригласить эксперта
Ваш ответ на вопрос

Войдите, чтобы написать ответ

Войти через центр авторизации
Похожие вопросы
Вакансии с Хабр Карьеры
Python developer
Bell Integrator
До 350 000 ₽
Python developer
Greenway Global • Новосибирск
от 150 000 ₽
TeamLead Python
AGIMA • Москва
До 350 000 ₽
Ещё вакансии
Заказы с Хабр Фриланса
Установка OS debian 11, raid 5, raid 0
04 мая 2024, в 22:32
2000 руб./за проект
Отправка сообщений в телеграм
04 мая 2024, в 22:17
12000 руб./за проект
Рассылка сообщений в телеграм
04 мая 2024, в 22:17
10000 руб./за проект
Ещё заказы