¿Cómo detener gdalwarp creación mundial salidas cerca de la fecha límite?

Question

Estoy usando gdalwarp para manipular SRTM azulejos cerca de la línea de cambio de fecha (es decir, 180°, también conocido como el antimeridian). SRTM azulejos tienen una muy leve (1/2 píxel) se superponen con el meridiano. Usted puede ver esta usando gdalinfo:

gdalinfo S16W180.hgt
Driver: SRTMHGT/SRTMHGT File Format
Files: S16W180.hgt
Size is 1201, 1201
[...]
Lower Left  (-180.0004167, -16.0004167) (180d 0' 1.50"W, 16d 0' 1.50"S)
Upper Right (-178.9995833, -14.9995833) (178d59'58.50"W, 14d59'58.50"S)
[...]

Para que la fuente se extiende la línea de cambio de fecha por una pequeña cantidad.

Esto causa problemas con gdalwarp, lo que termina generando un enorme globo-que abarca las salidas.

gdalwarp -t_srs "epsg:900913" S16W180.hgt test.tif
gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 1703, 5
[...]
Lower Left  (-20037508.330,-1806798.473) (180d 0' 0.00"W, 16d 7'13.00"S)
Upper Right (20032839.451,-1689152.120) (179d57'29.01"E, 15d 5'45.84"S)

Nota: las longitudes lapso de (casi) todo el mundo, y también el número de líneas es inesperadamente pequeño (5)

Es esto un error en gdalwarp? Si no, ¿cuáles son las opciones correctas para pasar a gdalwarp para obtener una salida sensata?

Answer 1

Funciona en dos pasos:

gdalwarp -te -180 -16 -179 -15 s16W180.hgt test.tif
gdalwarp -t_srs "epsg:3857" test.tif out.tif

El primer comando se inicia el píxel extra la mitad en el lado equivocado del meridiano 180°. Obtendrá un archivo de salida que es de 1178 P x 1222 L.

Por otra parte, con gdal_translate:

gdal_translate -a_ullr -180 -15 -179 -16 S16W180.hgt test2.tif
gdalwarp -t_srs "epsg:3857" test2.tif out2.tif

Crear un archivo de salida que es de 1179 P x 1223 L.

Answer 2

Una solución fácil sería para especificar el sistema de coordenadas "manualmente" como un PROJ cadena. Esto le permite utilizar el +over interruptor que desactiva la envoltura de la antimeridian:

gdalwarp -t_srs \
    "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0 \
        +over +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null \
        +wktext +lon_wrap=-180 +no_defs" \
    S16W180.hgt test.tif

Cuando hago eso y, a continuación, hacer gdalinfo sobre el resultado, me sale esto:

Corner Coordinates:
Upper Left  (-20037554.726,-1689152.120) (179d59'58.50"E, 14d59'58.50"S)
Lower Left  (-20037554.726,-1804766.925) (179d59'58.50"E, 16d 0' 1.37"S)
Upper Right (-19926099.407,-1689152.120) (178d59'57.11"W, 14d59'58.50"S)
Lower Right (-19926099.407,-1804766.925) (178d59'57.11"W, 16d 0' 1.37"S)
Center      (-19981827.066,-1746959.523) (179d29'59.30"W, 15d30' 2.12"S)

Tengo el PROJ cadena (sin +over) a partir de la observación de la salida original de gdalinfo. Fue incluido en un EXTENSION[...] de bloque del sistema de coordenadas.

Answer 3

Este es el problema en la biblioteca GDAL. Parece que GDALSuggestedWarpOutput() está dando salida extraña para la anchura y altura del archivo de salida.

No he encontrado una manera de trabajar alrededor de esto sin embargo.