Thursday, 2 November 2017

coordinate system - How to reproject raster from 0 360 to -180 180 with cutting 180 meridian


I have a geotiff raster image that has a coordinate system with longtitudes from 0 to 360. The horizontal center of the image is 180 longtitude. See image below:


enter image description here


I want to transform it to EPSG:4326 SRS with -180 180 longtitude range. And I want the center of the image to be at Greenwich meridian (0). I guess this srs is very widely used. I expect the result look like this:


enter image description here


So I use a gdalwarp command to reproject:



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

But I only get a tiff with bigger dimensions (more pixels) and EPSG:4326 metadata. The image itself looks the same, as the initial one. But I expect it to swap the hemispheres.


The question is - how do I gdalwarp an image to be strictly -180 180 EPSG:4326 with the center in 0 longtitude?


This is gdalinfo of my initial file:


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)

This is gdalinfo after 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)

Answer



You could explicitly set the output coordinate range using the target extent option to gdalwarp (ie. "-te -180 -90 180 90") but you can also use the CENTER_LONG configuration option to force rewrapping around a new central longitude. Something like this:


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

Note also the "SOURCE_EXTRA=1000" warp option. When doing rewrapping the source rectangle computation will get confused around the longitude interruption and lose some imagery. This option says pull in some extra. Without it you will see a data gap near the prime meridian.


PS. setting a prime meridian of 180dW as you did is not a good idea IMHO.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...