Bagaimana cara menghasilkan grid spasial dari raster?

9

Saya membutuhkan Spatial Grid sebagai master grid untuk beragam peta tematik. Bagaimana cara menghasilkan Spatial Grid dari raster yang membuang semua piksel NA?

Janbob Squarebrains
sumber
6
Beri kami beberapa memo di sini. Sedikit kode untuk membuat raster dan apa yang Anda harapkan darinya?
Spacedman

Jawaban:

14

Anda bisa mendapatkan semua koordinat sel non-NA dalam raster dengan:

r = raster(matrix(runif(20),5,4))
r[r>.5]=NA

coordinates(r)[!is.na(values(r)),]
          x   y
 [1,] 0.375 0.7
 [2,] 0.125 0.5
 [3,] 0.375 0.5
 [4,] 0.625 0.5
 [5,] 0.875 0.5
 [6,] 0.125 0.3
 [7,] 0.375 0.3
 [8,] 0.625 0.3
 [9,] 0.375 0.1
[10,] 0.875 0.1

mereka adalah sel-sel yang bukan NA. Anda mungkin dapat memberi ini ke SpatialPixels

SpatialPixels(SpatialPoints(coordinates(r)[!is.na(values(r)),]))
Object of class SpatialPixels
Grid topology:
  cellcentre.offset cellsize cells.dim
x             0.125     0.25         4
y             0.100     0.20         4
SpatialPoints:
          x   y
 [1,] 0.375 0.7
 [2,] 0.125 0.5
 [3,] 0.375 0.5
 [4,] 0.625 0.5
 [5,] 0.875 0.5
 [6,] 0.125 0.3
 [7,] 0.375 0.3
 [8,] 0.625 0.3
 [9,] 0.375 0.1
[10,] 0.875 0.1
Coordinate Reference System (CRS) arguments: NA 

Meskipun secara pribadi apa pun di grid saya tetap sebagai raster.

Saya tidak benar-benar yakin apa yang Anda inginkan - SpatialGridobjek didefinisikan grid persegi panjang penuh, jadi satu tanpa piksel NA tidak masuk akal.

Spacedman
sumber
terima kasih sobat, saya juga tidak yakin mau mau apa. Saya mencoba mengembangkan alur kerja untuk menghasilkan peta raster yang identik (dalam hal resolusi, berkoordinasi, crs, ...) dari berbagai data titik spasial. rasterize bukan pilihan karena hamburan titik spasial. Perlu menggunakan iwd atau serupa. Anda benar, kisi adalah segi empat - yang saya butuhkan sebagai master template mungkin adalah SpatialPointsDataFrame. yang dapat disambungkan. S.
Janbob Squarebrains
Anda mungkin menginginkan master raster (ha!) Dan kemudian membuat semua kisi-kisi berikutnya secara bersamaan.
Radar
2
Sepakat. Jika ada yang grid, raster hampir selalu jawabannya. Entah itu, atau mengarungi sekitar 20 halaman Analisis Data Spasial Terapan di R untuk mengetahui cara melakukan SpatialPixels. Saya memberikan pilihan itu kepada mahasiswa PhD hari ini. Dia memilih ... dengan bijak!
Spacedman
20

Untuk mengubah RasterLayer menjadi objek Spasial * (Kotak atau Piksel) Anda dapat menggunakan fungsi paksaan "sebagai"

library(raster)
r <- raster(matrix(runif(20),5,4))
r[r>.5] <- NA
g <- as(r, 'SpatialGridDataFrame')
p <- as(r, 'SpatialPixels')

plot(r)
points(p)
Robert Hijmans
sumber
7

Dua persyaratan Anda tampaknya tentang hal-hal yang berbeda:

1) Beberapa jenis templat raster grid yang andal.

2) Kisi jarang yang tidak secara eksplisit menyimpan sel yang hilang.

sp :: GridTopology menyediakan yang pertama, itu hanya deskripsi grid berdasarkan koordinat sel kiri bawah (cellcentre.offset), jarak sel (cellsize), dan dimensi grid (cells.dim).

Kelas sp :: SpatialPixelsDataFrame memungkinkan Anda untuk menyimpan grid jarang, tetapi dengan sendirinya menyimpan lebih dari "templat" - ia juga menyimpan setiap koordinat secara eksplisit. Ini karena ia melakukan dua pekerjaan, satu memungkinkan Anda untuk mempertahankan koordinat asli yang berasal dari grid dan mungkin sedikit tidak teratur, dua memungkinkan Anda menyimpan hanya sel-sel yang memiliki nilai yang valid. (Bisa dibilang * dua tujuan ini seharusnya dipisahkan, tapi itu cerita lain).

Saya tidak berpikir paket raster memiliki analog eksplisit dengan GridTopology, tetapi Anda bisa mendapatkan komponen untuk "roll your own":

library(raster)
r1 <- raster(nrows=108, ncols=21, xmn=0, xmx=10)

## "cellsize"
res(r1)
## [1] 0.4761905 1.6666667

## extreme cell corners (just a different convention to sp's cellcentre)
extent(r1)
class       : Extent 
xmin        : 0 
xmax        : 10 
ymin        : -90 
ymax        : 90 

## we can also use bbox to get the same thing
bbox(r1)
min max
s1   0  10
s2 -90  90

## grid dimensions (including number of attributes/layers as 3rd "dim")
dim(r1)
## [1] 108  21   1

Menghubungkan semua ini bersama-sama, kita bisa beralih dari raster ke sp:

GridTopology(bbox(r1)[,1] + res(r1)/2, res(r1), dim(r1)[2:1])

(Perhatikan bagaimana dimensi harus dibalik). Cara lain yang lebih sederhana adalah dengan memaksa SpatialGrid dan menggunakan sp's getGridTopology, meskipun ini lebih mahal karena koordinat akhirnya dihasilkan di sepanjang jalan:

getGridTopology(as(r1, "SpatialGrid"))

Ketiga bagian dari "topologi" raster tidak semuanya diperlukan, karena beberapa di antaranya berlebihan tetapi tidak ada cara lain untuk menangkap semuanya sebagai satu objek - kecuali, raster yang dibuat di atas adalah "kosong" sehingga dapat melakukan pekerjaan yang dilakukan GridTopology untuk sp. Saya tidak yakin tentang detail "kosong" itu, tetapi jelas tidak secara eksplisit menyimpan slot "data" dan lebih kecil daripada dengan nilai di dalamnya. Paket raster pada umumnya melakukan yang terbaik untuk menjaga penggunaan memori seminimal mungkin, dan karenanya Anda tidak perlu khawatir benar-benar "jarang".

Itu mungkin bisa membantu menjelaskan sedikit lagi, saya tahu saya tumpang tindih dengan jawaban Spacedman, tetapi masih belum jelas apa yang Anda maksud dalam pertanyaan itu.

  • (Dapat diperdebatkan, karena Anda dapat menyimpan sel jarang dengan hanya menyimpan indeks mereka daripada koordinat eksplisit, dan awalnya koordinat sel "sedikit tidak teratur" hanya dapat disimpan sebagai atribut jika Anda benar-benar menginginkannya. Baik sp atau raster tidak berurusan dengan raster yang tidak teratur sama sekali, bahkan bukan kasus bujursangkar sederhana - ini sesuai dengan sebagian besar alat SIG, tetapi sangat disayangkan karena ini cukup umum dalam format seperti NetCDF dan ditangani oleh grafis R :: :: fungsi gambar yang baik (meskipun tidak jika opsi rasterImage baru-baru ini digunakan) .)
mdsumner
sumber
Terima kasih atas pernyataan yang sangat instruktif, yang perlu saya cerna. Semoga memungkinkan saya untuk menyelesaikan masalah saya atau mengirim pertanyaan yang lebih eksplisit!
Janbob Squarebrains
Juga @Spacedman: Saya lebih tepatnya mendefinisikan latar belakang pertanyaan saya di
Janbob Squarebrains