Replicando el resultado de gdalwarp usando enlaces de Python gdal

17

Estoy intentando volver a proyectar / remuestrear con los enlaces de python de GDAL, pero obtengo resultados ligeramente diferentes en comparación con los de la utilidad de línea de comandos gdalwarp .

Consulte la actualización a continuación para ver un ejemplo más breve

Este script ilustra el enfoque de 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 con la salida de:

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

Se diferencian en tamaño (por 2 filas y 1 columna), así como con algunos valores de píxeles diferentes cerca de los bordes.

Vea la superposición transparente de test_epsg3573.tif y test_epsg3573_gdalwarp.tif a continuación. Si las imágenes fueran idénticas, solo habría píxeles blancos y negros, no gris.

ProbadoconPython2.7.8,GDAL1.11.1,Numpy1.9.1

Actualizar:

Aquíhayunejemplomuchomáscorto.Parecequeestonosedebealamejoradelmuestreoyaquelosiguientetambiénproduceresultadosinconsistentescongdalwarp

fromosgeoimportosr,gdalimportnumpy#Createsyntheticdatagtiff_drv=gdal.GetDriverByName('GTiff')w,h=512,512raster=numpy.zeros((w,h),dtype=numpy.uint8)raster[::w/10,:]=255raster[:,::h/10]=255top_left=(-109764,215677)pixel_size=45src_srs=osr.SpatialReference()src_srs.ImportFromEPSG(3413)src_geotran=[top_left[0],pixel_size,0,top_left[1],0,-pixel_size]rows,cols=raster.shapesrc_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)#ReprojecttoEPSG:3573dest_srs=osr.SpatialReference()dest_srs.ImportFromEPSG(3573)int_ds=gdal.AutoCreateWarpedVRT(src_ds,src_srs.ExportToWkt(),dest_srs.ExportToWkt())#Makedestdatasetdest_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())#Closedatasetssrc_ds=Nonedest_ds=None

Yestaeslallamadagdalwarpqueesperoquesealamisma,peronoloes:

gdalwarp-s_srsEPSG:3413-t_srsEPSG:3573-ofGTifftest_epsg3413.tiftest_epsg3573_gdalwarp.tif

Laimagendeabajomuestracadaimagenbinariaresultantesuperpuestaal50%detransparencia.Lospíxelesgrisesclarossoninconsistenciasentrelosdosresultados.

    
pregunta Bruce Wallin 23.03.2015 - 06:06

1 respuesta

13

Obtengo los mismos resultados que gdalwarp de gdal.AutoCreateWarpedVRT si establezco el umbral de error en 0.125 para que coincida con el valor predeterminado (-et) en gdalwarp . Alternativamente, puedes configurar -et 0.0 en tu llamada a gdalwarp para que coincida con el valor predeterminado en gdal.AutoCreateWarpedVRT .

Ejemplo

Crea una referencia para comparar con:

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

Ejecute la proyección en Python (basado en el código de la función "warp_27 () en GDAL autotest suite ):

# 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" 
    
respondido por el Luke 24.03.2015 - 00:11

Lea otras preguntas en las etiquetas