@zedxxx

Как работать с растрами в памяти в 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)

    ... ???
  • Вопрос задан
  • 3468 просмотров
Решения вопроса 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) отдельно, иначе получается каша...
Ответ написан
Комментировать
Пригласить эксперта
Ваш ответ на вопрос

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

Войти через центр авторизации
Похожие вопросы
Rhonda Software Владивосток
от 60 000 до 120 000 ₽
AKCP Новосибирск
от 180 000 ₽
Сбер Москва
от 170 000 ₽
28 нояб. 2021, в 21:01
3000 руб./за проект
28 нояб. 2021, в 20:54
1000 руб./за проект
28 нояб. 2021, в 19:47
500 руб./за проект