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

Question

[zed]

Стоит задача перепроецирования растров на лету (в памяти) средствами 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)

    ... ???

Answer 1

[zed]

Придумал загрузку и выгрузку растров через 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) отдельно, иначе получается каша...