Suite

Reprojection raster de 0 360 à -180 180 en coupant le méridien 180 à l'aide de gdalwarp

Reprojection raster de 0 360 à -180 180 en coupant le méridien 180 à l'aide de gdalwarp


J'ai une image raster geotiff qui a un système de coordonnées avec des longitudes de 0 à 360. Le centre horizontal de l'image est de 180 longitude. Voir image ci-dessous :

Je veux le transformer en EPSG:4326 SRS avec une plage de longitude -180 180. Et je veux que le centre de l'image soit au méridien de Greenwich (0). Je suppose que ce srs est très largement utilisé. Je m'attends à ce que le résultat ressemble à ceci:

J'utilise donc une commande gdalwarp pour reprojeter :

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

Mais je n'obtiens qu'un tiff avec des dimensions plus grandes (plus de pixels) et des métadonnées EPSG:4326. L'image elle-même est la même que l'image initiale. Mais je m'attends à ce qu'il échange les hémisphères.

Comment gdalwarp une image pour qu'elle soit strictement -180 180 EPSG:4326 avec le centre en longitude 0 ?

Voici gdalinfo de mon fichier initial :

Origine = (-0,102272598067084,89.946211604095552) Taille de pixel = (0.204545196134167,-0.204423208191126) Métadonnées : AREA_OR_POINT=Area Image Structure Métadonnées : INTERLEAVE=BAND Coordonnées d'angle : Supérieur gauche ( -0.1022726, 89.9462116) (0d 6' 8.18"W, 89d56' 46.36"N) Inférieur gauche (-0.1022726, -89.9462116) ( 0d 6' 8.18"W, 89d56'46.36"S) Supérieur droit ( 359.897, 89.946) (359d53'50.18"E, 89d56'46.36"N) Inférieur droit ( 359.897, -89.946) (359d53'50.18"E, 89d56'46.36"S) Centre ( 179.8975000, -0.0000000) (179d53'51.00"E, 0d 0' 0.00"S)

C'est gdalinfo après gdalwarp

Origine = (-180.102727401932952,89.946211604095552) Taille de pixel = (0.091397622896436,-0.091420837939082) Métadonnées : AREA_OR_POINT=Métadonnées de structure d'image de zone : INTERLEAVE=PIXEL Coordonnées de coin : en haut à gauche (-180.1027274, 89.9462116) (180d 6' 9.82'"W, 46.36"N) Inférieur gauche (-180.1027274, -89.9699975) (180d 6' 9.82"W, 89d58'11.99"S) Supérieur droit ( 179.821116, 89.9462116) (179d49'16.00"E, 89d56'46.36"N) Inférieur droit ( 179.821116, -89.9699975) (179d49'16.00"E, 89d58'11.99"S) Centre ( -0.1408079, -0.0118929) ( 0d 8'26.91"W, 0d 0'42.81"S)

Vous pouvez définir explicitement la plage de coordonnées de sortie à l'aide de l'option d'étendue cible sur gdalwarp (c'est-à-dire "-te -180 -90 180 90"), mais vous pouvez également utiliser l'option de configuration CENTER_LONG pour forcer le réenroulement autour d'une nouvelle longitude centrale. Quelque chose comme ça:

gdalwarp -s_srs "+proj=longlat +ellps=WGS84" -t_srs WGS84 ~/0_360.tif 180.tif -wo SOURCE_EXTRA=1000  --config CENTER_LONG 0

Notez également l'option de déformation "SOURCE_EXTRA=1000". Lors du réemballage, le calcul du rectangle source sera confus autour de l'interruption de la longitude et perdra certaines images. Cette option dit de tirer un peu plus. Sans cela, vous verrez une lacune de données près du méridien principal.

Je pense que définir un premier méridien de 180 dW comme vous l'avez fait n'est pas une bonne idée.


Fondamentalement, vous devez couper le raster en deux parties et les reconstituer avec un nouveau décalage/échelle.

Il y a un exemple ici de la façon de procéder de [-180,180] à [0,360] avec gdal_translate et le pilote VRT : http://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial

Scannez jusqu'au "Tutoriel de 5 minutes" et les détails se trouvent sous "Fichiers virtuels". Il devrait être assez simple de modifier l'exemple en conséquence.


Cela peut être fait dans R avec une ligne de code en utilisant letournerfonctionner avec lerasterpaquet.

chargez d'abord le package raster.

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

Si vous souhaitez simplement afficher le raster dans QGIS, vous pouvez définir une projection personnalisée avec le paramètre +lon_wrap=180.

Ma compréhension de ceci est que, par défaut, proj4 enveloppe les latitudes de 0 -> 360 à -180 -> 180. +lon_wrap=180 annulera effectivement cet enroulement et affichera les latitudes entre 180 et 360 dans l'hémisphère occidental.

L'option +over devrait désactiver complètement le wrapping, mais - du moins dans mon cas - le raster ne s'affichait pas correctement lorsque cette option était utilisée.

Voir http://proj4.org/parameters.html#lon-wrap-over-longitude-wrapping pour plus d'informations.


Voici une fonction que j'ai créée pour reprojeter un seul tableau de valeurs de grille à l'aide de JavaScript de 0-360 à -180-180.

let xstart = 180 / xres //xres est le nombre de valeurs pour 1 degré pour (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))) }