Cómo reproyectar ráster de 0 360 a -180 180 con corte de meridiano 180

31

Tengo una imagen ráster geotiff que tiene un sistema de coordenadas con longitudes de 0 a 360. El centro horizontal de la imagen es 180 de longitud. Ver imagen a continuación:

ingrese la descripción de la imagen aquí

Quiero transformarlo a EPSG: 4326 SRS con -180 180 rango de longitud. Y quiero que el centro de la imagen esté en el meridiano de Greenwich (0). Supongo que este srs es muy utilizado. Espero que el resultado se vea así:

ingrese la descripción de la imagen aquí

Entonces uso un comando gdalwarp para reproyectar:

gdalwarp -s_srs '+proj=latlong +datum=WGS84 +pm=180dW' -t_srs EPSG:4326 test_col.tif test_4326.tif

Pero solo obtengo un tiff con dimensiones más grandes (más píxeles) y EPSG: 4326 metadatos. La imagen en sí se ve igual que la inicial. Pero espero que cambie los hemisferios.

La pregunta es: ¿cómo puedo deformar una imagen para que sea estrictamente -180 180 EPSG: 4326 con el centro en longitud 0?

Este es gdalinfo de mi archivo inicial:

Origin = (-0.102272598067084,89.946211604095552)
Pixel Size = (0.204545196134167,-0.204423208191126)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -0.1022726,  89.9462116) (  0d 6' 8.18"W, 89d56'46.36"N)
Lower Left  (  -0.1022726, -89.9462116) (  0d 6' 8.18"W, 89d56'46.36"S)
Upper Right (     359.897,      89.946) (359d53'50.18"E, 89d56'46.36"N)
Lower Right (     359.897,     -89.946) (359d53'50.18"E, 89d56'46.36"S)
Center      ( 179.8975000,  -0.0000000) (179d53'51.00"E,  0d 0' 0.00"S)

Esto es gdalinfo después de gdalwarp

Origin = (-180.102727401932952,89.946211604095552)
Pixel Size = (0.091397622896436,-0.091420837939082)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.1027274,  89.9462116) (180d 6' 9.82"W, 89d56'46.36"N)
Lower Left  (-180.1027274, -89.9699975) (180d 6' 9.82"W, 89d58'11.99"S)
Upper Right ( 179.8211116,  89.9462116) (179d49'16.00"E, 89d56'46.36"N)
Lower Right ( 179.8211116, -89.9699975) (179d49'16.00"E, 89d58'11.99"S)
Center      (  -0.1408079,  -0.0118929) (  0d 8'26.91"W,  0d 0'42.81"S)
nextstopsun
fuente
Acerca de la resolución diferente, ¿has intentado agregar la -tr xres yresbandera?
nickves

Respuestas:

21

Puede establecer explícitamente el rango de coordenadas de salida utilizando la opción de extensión de destino para gdalwarp (es decir, "-te -180 -90 180 90") pero también puede usar la opción de configuración CENTER_LONG para forzar el reenvolvimiento alrededor de una nueva longitud central. Algo como esto:

  gdalwarp -t_srs WGS84 ~/0_360.tif 180.tif  -wo SOURCE_EXTRA=1000 \
           --config CENTER_LONG 0

Tenga en cuenta también la opción de deformación "SOURCE_EXTRA = 1000". Al volver a envolver el cálculo del rectángulo de origen se confundirá con la interrupción de longitud y perderá algunas imágenes. Esta opción dice tirar un poco más. Sin él, verá una brecha de datos cerca del meridiano principal.

PD. establecer un meridiano principal de 180dW como lo hizo no es una buena idea en mi humilde opinión.

Frank Warmerdam
fuente
1
hmm, --config CENTER_LONG 0no hace nada, el resultado es el mismo ráster. ¿Algo que extraño aquí? Ejecutando en GDAL versión 2.2.3.
jurajb
6

Básicamente, debe cortar el ráster en dos partes y juntarlas nuevamente con un nuevo desplazamiento / escala.

Aquí hay un ejemplo de cómo hacerlo desde [-180,180] a [0,360] con gdal_translate y el controlador VRT: http://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial

Escanee hacia abajo hasta el "tutorial de 5 minutos" y los detalles se encuentran en "Archivos virtuales". Debería ser lo suficientemente simple como para modificar el ejemplo a la medida.

mdsumner
fuente
2

Esto se puede hacer en R con una línea de código usando la rotatefunción con el rasterpaquete.

library(raster)
your_raster <- raster("path/to/raster.tif")
rotated_raster <- rotate(your_raster)
SoilSciGuy
fuente
1

Si solo desea ver el ráster en QGIS, puede establecer una proyección personalizada con el parámetro + lon_wrap = 180.

Tengo entendido que, por defecto, proj4 ajusta las latitudes de 0 -> 360 a -180 -> 180. + lon_wrap = 180 cancelará efectivamente este ajuste y mostrará latitudes entre 180 y 360 en el hemisferio occidental.

La opción + over debería deshabilitar la envoltura por completo, pero, al menos en mi caso, el ráster no se mostró correctamente cuando se usó esa opción.

Consulte http://proj4.org/parameters.html#lon-wrap-over-longitude-wrapping para obtener más información.


fuente
0

Aquí hay una función que construí para reproyectar una única matriz tenue de valores de cuadrícula usando javascript de 0-360 a -180-180. Espero que pueda ser de ayuda para alguien.

  let xstart = 180 / xres //xres is the number of values per 1 degree
  for (let y = 0; y < data.height; y++) {
    let index = (y * data.width) + 1,
    start = index + xstart,
    end = index + data.width
    array.splice(index, 0, ...array.splice(start, (end - start)))
  }
maeneak
fuente