Cara memproyeksikan ulang raster dari 0 360 ke -180 180 dengan memotong 180 meridian

31

Saya memiliki gambar raster geotiff yang memiliki sistem koordinat dengan longtitude dari 0 hingga 360. Pusat horizontal gambar adalah 180 longtitude. Lihat gambar di bawah ini:

masukkan deskripsi gambar di sini

Saya ingin mengubahnya menjadi EPSG: 4326 SRS dengan -180 180 rentang longtitude. Dan saya ingin bagian tengah gambar berada di Greenwich meridian (0). Saya kira srs ini sangat banyak digunakan. Saya berharap hasilnya terlihat seperti ini:

masukkan deskripsi gambar di sini

Jadi saya menggunakan perintah gdalwarp untuk memproyeksi ulang:

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

Tapi saya hanya mendapatkan tiff dengan dimensi lebih besar (lebih banyak piksel) dan EPSG: 4326 metadata. Gambar itu sendiri terlihat sama, seperti yang pertama. Tapi saya berharap itu untuk menukar belahan otak.

Pertanyaannya adalah - bagaimana cara membuat gambar menjadi ketat -180 180 EPSG: 4326 dengan pusat dalam 0 longtitude?

Ini adalah gdalinfo dari file awal saya:

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)

Ini adalah gdalinfo setelah 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
sumber
Tentang resolusi yang berbeda, sudahkah Anda mencoba menambahkan -tr xres yresbenderanya?
nickves

Jawaban:

21

Anda dapat secara eksplisit mengatur rentang koordinat keluaran menggunakan opsi batas target ke gdalwarp (mis. "-Te -180 -90 180 90") tetapi Anda juga dapat menggunakan opsi konfigurasi CENTER_LONG untuk memaksa membungkus kembali sekitar bujur pusat baru. Sesuatu seperti ini:

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

Perhatikan juga opsi warp "SOURCE_EXTRA = 1000". Saat melakukan rewrapping, perhitungan persegi panjang sumber akan kebingungan di sekitar gangguan garis bujur dan kehilangan beberapa citra. Opsi ini mengatakan menarik beberapa tambahan. Tanpanya Anda akan melihat celah data di dekat meridian utama.

PS. menetapkan meridian utama 180dW seperti yang Anda lakukan bukanlah IMHO ide yang baik.

Frank Warmerdam
sumber
1
hmm, --config CENTER_LONG 0tidak apa-apa, hasilnya adalah raster yang sama. Ada yang saya lewatkan di sini? Berjalan di GDAL versi 2.2.3.
jurajb
6

Pada dasarnya Anda perlu memotong raster menjadi dua bagian dan menyatukannya kembali dengan offset / skala baru.

Ada contoh di sini tentang bagaimana melakukan itu dari [-180.180] ke [0,360] dengan gdal_translate dan driver VRT: http://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial

Pindai ke "tutorial 5 menit" dan detailnya ada di bawah "File Virtual". Itu harus cukup sederhana untuk memodifikasi contoh yang sesuai.

mdsumner
sumber
2

Ini dapat dilakukan dalam R dengan satu baris kode menggunakan rotatefungsi dengan rasterpaket.

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

Jika Anda hanya ingin melihat raster di QGIS, Anda dapat mengatur proyeksi kustom dengan parameter + lon_wrap = 180.

Pemahaman saya tentang hal ini adalah, secara default, proj4 membungkus garis lintang dari 0 -> 360 hingga -180 -> 180. + lon_wrap = 180 secara efektif akan membatalkan pembungkus ini, dan menampilkan garis lintang antara 180 dan 360 di belahan bumi Barat.

Opsi + lebih harus menonaktifkan pembungkus sama sekali, tetapi - setidaknya dalam kasus saya - raster tidak ditampilkan dengan benar ketika opsi itu digunakan.

Lihat http://proj4.org/parameters.html#lon-wrap-over-longitude-wrapping untuk informasi lebih lanjut.


sumber
0

Berikut adalah fungsi yang saya buat untuk memproyeksikan ulang array redup tunggal nilai grid menggunakan javascript dari 0-360 ke -180-180. Saya harap ini bisa membantu seseorang.

  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
sumber