Я пытаюсь перепроектировать / повторно сэмплировать с привязками Python GDAL, но получаю немного другие результаты по сравнению с результатами из утилиты командной строки gdalwarp
.
Смотрите обновление ниже для более короткого примера
Этот скрипт иллюстрирует подход 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
Сравните с выводом:
gdalwarp -s_srs EPSG:3413 -t_srs EPSG:3573 -tr 7 7 -r near -of GTiff test_epsg3413.tif test_epsg3573_gdalwarp.tif
Они различаются по размеру (на 2 строки и 1 столбец), а также с некоторыми различными значениями пикселей вблизи краев.
Смотрите прозрачный оверлей test_epsg3573.tif и test_epsg3573_gdalwarp.tif ниже. Если бы изображения были идентичными, то были бы только черные и белые пиксели, а не серые.
Протестировано с Python 2.7.8, GDAL 1.11.1, Numpy 1.9.1
Обновление :
Вот гораздо более короткий пример. Кажется, это не вызвано повышением частоты дискретизации, поскольку следующее также дает результаты, несовместимые сgdalwarp
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
И это вызов gdalwarp, который, я ожидаю, будет таким же, но это не так:
gdalwarp -s_srs EPSG:3413 -t_srs EPSG:3573 -of GTiff test_epsg3413.tif test_epsg3573_gdalwarp.tif
Изображение ниже показывает каждое полученное двоичное изображение, наложенное с прозрачностью 50%. Светло-серые пиксели - это несоответствия между двумя результатами.
gdal.AutoCreateWarpedVRT(source_file, source_srs_wkt, dest_srs_wkt)
?Ответы:
Я получаю те же результаты,
gdalwarp
отgdal.AutoCreateWarpedVRT
если установить пороговое значение ошибки 0,125 , чтобы соответствовать по умолчанию (-et) в gdalwarp . Кроме того, вы можете установить-et 0.0
в своем вызове,gdalwarp
чтобы соответствовать по умолчанию вgdal.AutoCreateWarpedVRT
.пример
Создать ссылку для сравнения с:
Запустите проекцию в Python (на основе кода из функции "warp_27 () в наборе автотестов GDAL ):
источник