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 ).
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.
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?Jawaban:
Dibutuhkan beberapa langkah:
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.
b. Konversi bingkai data ke SpatialPointsDataFrame menggunakan paket sp dan sesuatu seperti:
c. Konversikan ke sistem km reguler Anda dengan terlebih dahulu memberi tahu apa itu CRS, dan kemudian spTransformasi ke tujuan.
d. Beri tahu R bahwa ini dikisi:
Pada titik ini Anda akan mendapatkan kesalahan jika koordinat Anda tidak terletak pada kisi biasa yang bagus.
Sekarang gunakan paket raster untuk mengkonversi ke raster dan atur CRS-nya:
Sekarang lihat:
Sekarang tulis sebagai file geoTIFF menggunakan paket raster:
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 ...
sumber
r = Import[ "https://ebi-forecast.igb.illinois.edu/bety/miscanthusyield.csv", "Data"];
. Anda bisa mendapatkan plot cepat poin setelahnya melaluidata = Rest[r]; ListPlot[data[[;; , {3, 2}]]]
(atauListPointPlot3D[data[[;; , {3, 2, 4}]]]
). Untuk proyeksi ulang, mulailah dengan bantuanGeoGridPosition
, 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.Sejak pertanyaan terakhir dijawab, ada solusi yang jauh lebih mudah dengan menggunakan fungsi paket raster
rasterFromXYZ
yang merangkum semua langkah yang diperlukan (termasuk spesifikasi string CRS).sumber
/tmp/
itu ~ 19GB ketika saya kehabisan ruang disk.merge()
hasilnya bersama-sama.