Replicando o resultado do gdalwarp usando as ligações gdal Python

20

Estou tentando reprojetar / reamostrar com as ligações python GDAL, mas estou obtendo resultados ligeiramente diferentes em comparação aos do utilitário de linha de comando gdalwarp.

Veja a atualização abaixo para obter um exemplo mais curto

Este script ilustra a abordagem do Python:

from osgeo import osr, gdal
import numpy


def reproject_point(point, srs, target_srs):
    '''
    Reproject a pair of coordinates from one spatial reference system to
    another.
    '''
    transform = osr.CoordinateTransformation(srs, target_srs)
    (x, y, z) = transform.TransformPoint(*point)

    return (x, y)


def reproject_bbox(top_left, bottom_right, srs, dest_srs):
    x_min, y_max = top_left
    x_max, y_min = bottom_right
    corners = [
        (x_min, y_max),
        (x_max, y_max),
        (x_max, y_min),
        (x_min, y_min)]
    projected_corners = [reproject_point(crnr, srs, dest_srs)
                         for crnr in corners]

    dest_top_left = (min([crnr[0] for crnr in projected_corners]),
                     max([crnr[1] for crnr in projected_corners]))
    dest_bottom_right = (max([crnr[0] for crnr in projected_corners]),
                         min([crnr[1] for crnr in projected_corners]))

    return dest_top_left, dest_bottom_right


################################################################################
# Create synthetic data
gtiff_drv = gdal.GetDriverByName('GTiff')
w, h = 512, 512
raster = numpy.zeros((w, h), dtype=numpy.uint8)
raster[::w / 10, :] = 255
raster[:, ::h / 10] = 255
top_left = (-109764, 215677)
pixel_size = 45

src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(3413)

src_geotran = [top_left[0], pixel_size, 0,
               top_left[1], 0, -pixel_size]

rows, cols = raster.shape
src_ds = gtiff_drv.Create(
    'test_epsg3413.tif',
    cols, rows, 1,
    gdal.GDT_Byte)
src_ds.SetGeoTransform(src_geotran)
src_ds.SetProjection(src_srs.ExportToWkt())
src_ds.GetRasterBand(1).WriteArray(raster)


################################################################################
# Reproject to EPSG: 3573 and upsample to 7m
dest_pixel_size = 7

dest_srs = osr.SpatialReference()
dest_srs.ImportFromEPSG(3573)

# Calculate new bounds by re-projecting old corners
x_min, y_max = top_left
bottom_right = (x_min + cols * pixel_size,
                y_max - rows * pixel_size)
dest_top_left, dest_bottom_right = reproject_bbox(
    top_left, bottom_right,
    src_srs, dest_srs)

# Make dest dataset
x_min, y_max = dest_top_left
x_max, y_min = dest_bottom_right
new_rows = int((x_max - x_min) / float(dest_pixel_size))
new_cols = int((y_max - y_min) / float(dest_pixel_size))
dest_ds = gtiff_drv.Create(
    'test_epsg3573.tif',
    new_rows, new_cols, 1,
    gdal.GDT_Byte)
dest_geotran = (dest_top_left[0], dest_pixel_size, 0,
                dest_top_left[1], 0, -dest_pixel_size)
dest_ds.SetGeoTransform(dest_geotran)
dest_ds.SetProjection(dest_srs.ExportToWkt())

# Perform the projection/resampling
gdal.ReprojectImage(
    src_ds, dest_ds,
    src_srs.ExportToWkt(), dest_srs.ExportToWkt(),
    gdal.GRA_NearestNeighbour)

dest_data = dest_ds.GetRasterBand(1).ReadAsArray()

# Close datasets
src_ds = None
dest_ds = None

Compare com a saída de:

gdalwarp -s_srs EPSG:3413 -t_srs EPSG:3573 -tr 7 7 -r near -of GTiff test_epsg3413.tif test_epsg3573_gdalwarp.tif

Eles diferem em tamanho (por 2 linhas e 1 coluna), bem como por alguns valores de pixel diferentes nas bordas.

Consulte sobreposição transparente de test_epsg3573.tif e test_epsg3573_gdalwarp.tif abaixo. Se as imagens fossem idênticas, haveria apenas pixels preto e branco, sem cinza.

Sobreposição QGIS de test_epsg3573.tif e test_epsg3573_gdalwarp.tif

Testado com Python 2.7.8, GDAL 1.11.1, Numpy 1.9.1

Atualização :

Aqui está um exemplo muito mais curto. Isso parece não ser causado por upsampling, pois o seguinte também produz resultados inconsistentes comgdalwarp

from osgeo import osr, gdal
import numpy


# Create synthetic data
gtiff_drv = gdal.GetDriverByName('GTiff')
w, h = 512, 512
raster = numpy.zeros((w, h), dtype=numpy.uint8)
raster[::w / 10, :] = 255
raster[:, ::h / 10] = 255
top_left = (-109764, 215677)
pixel_size = 45

src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(3413)

src_geotran = [top_left[0], pixel_size, 0,
               top_left[1], 0, -pixel_size]

rows, cols = raster.shape
src_ds = gtiff_drv.Create(
    'test_epsg3413.tif',
    cols, rows, 1,
    gdal.GDT_Byte)
src_ds.SetGeoTransform(src_geotran)
src_ds.SetProjection(src_srs.ExportToWkt())
src_ds.GetRasterBand(1).WriteArray(raster)

# Reproject to EPSG: 3573
dest_srs = osr.SpatialReference()
dest_srs.ImportFromEPSG(3573)

int_ds = gdal.AutoCreateWarpedVRT(src_ds, src_srs.ExportToWkt(), dest_srs.ExportToWkt())

# Make dest dataset
dest_ds = gtiff_drv.Create(
    'test_epsg3573_avrt.tif',
    int_ds.RasterXSize, int_ds.RasterYSize, 1,
    gdal.GDT_Byte)
dest_ds.SetGeoTransform(int_ds.GetGeoTransform())
dest_ds.SetProjection(int_ds.GetProjection())
dest_ds.GetRasterBand(1).WriteArray(int_ds.GetRasterBand(1).ReadAsArray())

# Close datasets
src_ds = None
dest_ds = None

E este é o chamado de gdalwarp que espero ser o mesmo, mas não é:

gdalwarp -s_srs EPSG:3413 -t_srs EPSG:3573 -of GTiff test_epsg3413.tif test_epsg3573_gdalwarp.tif

A imagem abaixo mostra cada imagem binária resultante sobreposta com 50% de transparência. Os pixels cinza claro são inconsistências entre os dois resultados.

Inconsistência ilustrada no QGIS

python gdal gdalwarp numpy Bruce Wallin
fonte
1
Você já tentou gdal.AutoCreateWarpedVRT(source_file, source_srs_wkt, dest_srs_wkt)?
user2856
Obrigado Luke, não conhecia essa função. Tentei agora, mas alguns pixels ainda são diferentes entre os dois. Ou seja, as geotransformações e formas dos rasters são idênticas (quando não são amostradas), mas alguns pixels parecem ter uma nova amostra diferente. Isso pelo menos demonstra que o problema ainda está presente, mesmo quando não há amostragem.
23615 Bruce Wallin

Respostas:

16

Eu obtenho os mesmos resultados gdalwarpde gdal.AutoCreateWarpedVRTse eu definir o limite de erro como 0,125 para corresponder ao padrão (-et) no gdalwarp . Como alternativa, você pode definir -et 0.0na sua chamada como gdalwarpo padrão em gdal.AutoCreateWarpedVRT.

Exemplo

Crie uma referência para comparar:

gdalwarp -t_srs EPSG:4326 byte.tif warp_ref.tif

Execute a projeção em Python (com base no código da função "warp_27 () no conjunto de autoteste GDAL )):

# Open source dataset
src_ds = gdal.Open('byte.tif')

# Define target SRS
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(4326)
dst_wkt = dst_srs.ExportToWkt()

error_threshold = 0.125  # error threshold --> use same value as in gdalwarp
resampling = gdal.GRA_NearestNeighbour

# Call AutoCreateWarpedVRT() to fetch default values for target raster dimensions and geotransform
tmp_ds = gdal.AutoCreateWarpedVRT( src_ds,
                                   None, # src_wkt : left to default value --> will use the one from source
                                   dst_wkt,
                                   resampling,
                                   error_threshold )

# Create the final warped raster
dst_ds = gdal.GetDriverByName('GTiff').CreateCopy('warp_test.tif', tmp_ds)
dst_ds = None

# Check that we have the same result as produced by 'gdalwarp -rb -t_srs EPSG:4326 ....'

ref_ds = gdal.Open('warp_ref.tif')
ref_cs = ref_ds.GetRasterBand(1).Checksum()

ds = gdal.Open('warp_test.tif')
cs = ds1.GetRasterBand(1).Checksum()

if cs == ref_cs:
    print 'success, they match'
else:
    print "fail, they don't match"
user2856
fonte