Bagaimana saya bisa mengkonversi data dalam bentuk nilai lat, lon, menjadi file raster menggunakan R?

40

Saya memiliki sekumpulan nilai pada kisi km di benua AS. Kolomnya adalah "lintang", "bujur", dan "pengamatan", misalnya:

"lat"    "lon"     "yield"
 25.567  -120.347  3.6 
 25.832  -120.400  2.6
 26.097  -120.454  3.4
 26.363  -120.508  3.1
 26.630  -120.562  4.4

atau, sebagai bingkai data R:

mydata <- structure(list(lat = c(25.567, 25.832, 26.097, 26.363, 26.63), 
lon = c(-120.347, -120.4, -120.454, -120.508, -120.562), 
yield = c(3.6, 2.6, 3.4, 3.1, 4.4)), .Names = c("lat", 
"lon", "yield"), class = "data.frame", row.names = c(NA, -5L))

(set data lengkap dapat diunduh sebagai csv di sini )

Data tersebut adalah output dari model krop (dimaksudkan untuk berada di) kisi 30 km x 30 km (dari Miguez et al 2012 ).

masukkan deskripsi gambar di sini

Bagaimana saya bisa mengonversinya menjadi file raster dengan metadata terkait GIS seperti proyeksi peta?

Idealnya file tersebut akan menjadi file teks (ASCII?) Karena saya ingin itu menjadi platform dan perangkat lunak independen.

Abe
sumber
Sebagai CSV, ini sudah merupakan "file teks" di ASCII. Juga, karena tidak menggunakan proyeksi sama sekali, mungkin ada sedikit metadata yang relevan untuk ditambahkan (datum, kebanyakan). Bisakah Anda sedikit lebih spesifik tentang apa jenis output yang Anda cari dan apa yang ingin Anda lakukan dengannya?
whuber
Saya ingin membuatnya semudah mungkin bagi seseorang untuk menggunakan data dengan berbagai perangkat lunak pemetaan (ArcGIS, Google Maps, Grass, R, dll.) Sehingga memudahkan penggunaan kembali, misalnya dengan tidak memerlukan langkah-langkah konversi tambahan. Berdasarkan halaman Wikipedia dalam format file GIS, saya menyimpulkan 1) file "raster" harus memiliki nama belakang dengan garis lintang dan nama kolom garis bujur, seperti gambar dan bahwa 2) metadata harus mencakup informasi geografis (lokasi sudut, area yang tertutup oleh data).
Abe
Ini adalah salah satu referensi terbaik yang saya temui pada R dan GIS. Terima kasih banyak! Bisakah Anda memberikan csv lain lat dan panjang dengan proj4string yang benar? Saya akan sangat menghargai itu.
@Nandini Tidak yakin apa proj4string benar adalah, saya menduga lambert konformal:proj +proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +lon_0=-145.5 +lat_0=1.0 . Saya tidak yakin apa yang Anda minta sehubungan dengan file csv lain - bagaimana hal itu berbeda dari yang ditautkan ke dalam pertanyaan, atau yang akan dihasilkan oleh jawaban yang diterima?
Abe
bagi saya tidak bekerja! Saya tidak tahu harus mengenakan koordinat "x" dan "y" ke "(pts) = ~ x + y"

Jawaban:

44

Dibutuhkan beberapa langkah:

  1. Anda bilang itu grid 1km biasa, tapi itu artinya lat-long tidak teratur. Pertama, Anda perlu mengubahnya menjadi sistem koordinat grid biasa sehingga nilai X dan Y ditempatkan secara teratur.

    Sebuah. Bacakan ke R sebagai bingkai data, dengan kolom x, y, hasil.

    pts = read.table("file.csv",......)

    b. Konversi bingkai data ke SpatialPointsDataFrame menggunakan paket sp dan sesuatu seperti:

    library(sp)
    library(rgdal)
    coordinates(pts)=~x+y

    c. Konversikan ke sistem km reguler Anda dengan terlebih dahulu memberi tahu apa itu CRS, dan kemudian spTransformasi ke tujuan.

    proj4string(pts)=CRS("+init=epsg:4326") # set it to lat-long
    pts = spTransform(pts,CRS("insert your proj4 string here"))

    d. Beri tahu R bahwa ini dikisi:

    gridded(pts) = TRUE

    Pada titik ini Anda akan mendapatkan kesalahan jika koordinat Anda tidak terletak pada kisi biasa yang bagus.

  2. Sekarang gunakan paket raster untuk mengkonversi ke raster dan atur CRS-nya:

    r = raster(pts)
    projection(r) = CRS("insert your proj4 string here")
  3. Sekarang lihat:

    plot(r)
  4. Sekarang tulis sebagai file geoTIFF menggunakan paket raster:

    writeRaster(r,"pts.tif")

GeoTIFF ini harus dapat dibaca di semua paket GIS utama. Bagian yang hilang jelas di sini adalah string proj4 untuk mengkonversi: ini mungkin akan menjadi semacam sistem referensi UTM. Sulit untuk mengatakan tanpa lebih banyak data ...

Spacedman
sumber
+1 Terima kasih telah menyiapkan alur kerja. Perhatikan bahwa data tersedia di tautan yang disediakan dalam pertanyaan: lihat. Anda akan menemukan, sayangnya, bahwa beberapa asumsi Anda tentang mereka salah. (Secara khusus, saya diburu untuk setiap dokumentasi tentang proyeksi yang digunakan untuk membuat grid tetapi tidak menemukannya Dan itu adalah proyeksi yang aneh, seperti yang Anda lihat dengan memplot poin..)
whuber
Sangat dekat untuk menjadi sistem UTM, tetapi tidak ada yang saya coba cukup dekat dengan grid biasa untuk R untuk grid mereka. Saya setengah tergoda untuk mengulangi seluruh basis data epsg R ....
Spacedman
Itu akan menjadi tour de force nyata jika Anda bisa menemukan proyeksi seperti itu! Kuncinya adalah menemukan kriteria yang efektif dan efisien untuk menentukan kapan 7.000 poin ini cukup dekat untuk diletakkan di grid biasa (karena mungkin mereka mungkin tidak membentuk grid sempurna dalam proyeksi standar sama sekali). Untuk menjalankan cepat melalui database, itu harus cukup untuk membandingkan sejumlah kecil jarak, seperti jarak timur-barat di bagian utara grid ke jarak timur-barat di bagian selatan. Itu seharusnya menghilangkan sebagian besar kandidat dengan cepat.
whuber
3
Saya berlari melalui semua proyeksi (default) yang didukung oleh Mathematica 8. Ia menemukan proyeksi di mana poin-poinnya benar-benar jatuh pada grid: Alaska State Plane (1983) Zone 10! Ini adalah proyeksi Lambert Conformal Conic. Saya percaya ini adalah EPSG 26940 . Jika Anda memodifikasi ini untuk memusatkannya kira-kira di bujur -106, titik-titik membentuk kisi yang cukup bagus.
Whuber
1
Abe, maksud Anda membaca halaman Web? Itu tadi r = Import[ "https://ebi-forecast.igb.illinois.edu/bety/miscanthusyield.csv", "Data"];. Anda bisa mendapatkan plot cepat poin setelahnya melalui data = Rest[r]; ListPlot[data[[;; , {3, 2}]]](atau ListPointPlot3D[data[[;; , {3, 2, 4}]]]). Untuk proyeksi ulang, mulailah dengan bantuan GeoGridPosition, kemudian buat beberapa tebakan cerdas dan referensi silang untuk mencari tahu apa yang terjadi :-). BTW, penjelasan @ Spacedman benar-benar relevan: distorsi metrik dari 25 menjadi 49 derajat sama dengan cos (25) / cos (49) = 1,38; itu substansial.
whuber
29

Sejak pertanyaan terakhir dijawab, ada solusi yang jauh lebih mudah dengan menggunakan fungsi paket raster rasterFromXYZyang merangkum semua langkah yang diperlukan (termasuk spesifikasi string CRS).

library(raster)
rasterFromXYZ(mydata)
Lucas Fortini
sumber
1
Permintaan maaf kepada @Spacedman yang tak kenal lelah, yang telah sering membantu saya, tetapi saya pikir jawaban ini layak untuk mewarisi tanda centang yang ramah hijau.
geotheory
@geotheory Saya akan memilih jawaban ini, ini adalah fungsi yang hebat, tetapi tampaknya sangat lambat pada dataset yang saya gunakan (ditautkan ke dalam op)
Abe
1
... sebenarnya itu tersedak karena mengambil ~ 400KB file saya dan membuat file /tmp/itu ~ 19GB ketika saya kehabisan ruang disk.
Abe
Mungkin ada proses n-kuadrat di sana di suatu tempat. Anda mungkin dapat mengelompokkan data poin dengan kisi yang luas, rasterise masing-masing kelompok secara individu dan kemudian merge()hasilnya bersama-sama.
geoteori
Dengan segala hormat, tapi jawaban ini jauh lebih baik daripada Spacedman.
Hantu