Konversikan garis shapefile menjadi raster, nilai = total panjang garis dalam sel

14

Saya memiliki garis shapefile yang mewakili jaringan jalan. Saya ingin merasterisasi data ini, dengan nilai yang dihasilkan di raster menunjukkan panjang total garis yang termasuk dalam sel raster.

Data dalam proyeksi British National Grid sehingga unit akan menjadi meter.

Idealnya saya ingin melakukan operasi ini menggunakan R, dan saya menduga bahwa rasterizefungsi dari rasterpaket akan berperan dalam mencapai ini, saya hanya tidak bisa mengetahui apa fungsi yang seharusnya diterapkan.

JPD
sumber
Mungkin vignette('over', package = 'sp')bisa membantu.

Jawaban:

12

Mengikuti pertanyaan terakhir , Anda mungkin ingin menggunakan fungsionalitas yang ditawarkan oleh paket rgeos untuk menyelesaikan masalah Anda. Untuk alasan reproduktifitas, saya mengunduh bentuk jalan Tanzania dari DIVA-GIS dan memasukkannya ke direktori kerja saya saat ini. Untuk tugas yang akan datang, Anda akan membutuhkan tiga paket:

  • rgdal untuk penanganan data spasial umum
  • raster untuk rasterization dari data shapefile
  • rgeos untuk memeriksa persimpangan jalan dengan templat raster dan menghitung panjang jalan

Akibatnya, baris pertama Anda dapat terlihat seperti ini:

library(rgdal)
library(raster)
library(rgeos)

Setelah itu, Anda perlu mengimpor data shapefile. Perhatikan bahwa DIVA-GIS shapefile didistribusikan dalam EPSG: 4326, jadi saya akan memproyeksikan shapefile ke EPSG: 21037 (UTM 37S) untuk menangani meter daripada derajat.

roads <- readOGR(dsn = ".", layer = "TZA_roads")
roads_utm <- spTransform(roads, CRS("+init=epsg:21037"))

Untuk rasterisasi selanjutnya, Anda akan memerlukan templat raster yang mencakup batas spasial dari shapefile Anda. Templat raster terdiri dari 10 baris dan 10 kolom secara default, sehingga menghindari waktu komputasi yang terlalu luas.

roads_utm_rst <- raster(extent(roads_utm), crs = projection(roads_utm))

Sekarang templat sudah diatur, loop melalui semua sel raster (yang saat ini hanya terdiri dari nilai-nilai NA). Dengan menetapkan nilai '1' ke sel saat ini dan kemudian mengeksekusi rasterToPolygons, shapefile yang dihasilkan 'tmp_shp' secara otomatis menampung luasnya piksel yang saat ini diproses. gIntersectsmendeteksi apakah tingkat ini tumpang tindih dengan jalan. Jika tidak, fungsi akan mengembalikan nilai '0'. Jika tidak, bentuk jalan dipotong oleh sel saat ini dan total panjang 'SpatialLines' dalam sel itu sedang dihitung menggunakan gLength.

lengths <- sapply(1:ncell(roads_utm_rst), function(i) {
  tmp_rst <- roads_utm_rst
  tmp_rst[i] <- 1
  tmp_shp <- rasterToPolygons(tmp_rst)

  if (gIntersects(roads_utm, tmp_shp)) {
    roads_utm_crp <- crop(roads_utm, tmp_shp)
    roads_utm_crp_length <- gLength(roads_utm_crp)
    return(roads_utm_crp_length)
  } else {
    return(0)
  }
})

Terakhir, Anda dapat memasukkan panjang yang dihitung (yang dikonversi menjadi kilometer) ke dalam templat raster dan secara visual memverifikasi hasil Anda.

roads_utm_rst[] <- lengths / 1000

library(RColorBrewer)
spplot(roads_utm_rst, scales = list(draw = TRUE), xlab = "x", ylab = "y", 
       col.regions = colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout = list("sp.lines", roads_utm), 
       par.settings = list(fontsize = list(text = 15)), at = seq(0, 1800, 200))

lines2raster

fdetsch
sumber
Ini luar biasa! Saya mengubah sapply()ke pbsapply()dan digunakan argumen klaster cl = detectCores()-1. Sekarang saya dapat menjalankan contoh ini secara paralel!
philiporlando
11

Di bawah ini dimodifikasi dari solusi Jeffrey Evans. Solusi ini jauh lebih cepat karena tidak menggunakan rasterize

library(raster)
library(rgdal)
library(rgeos)

roads <- shapefile("TZA_roads.shp")
roads <- spTransform(roads, CRS("+proj=utm +zone=37 +south +datum=WGS84"))
rs <- raster(extent(roads), crs=projection(roads))
rs[] <- 1:ncell(rs)

# Intersect lines with raster "polygons" and add length to new lines segments
rsp <- rasterToPolygons(rs)
rp <- intersect(roads, rsp)
rp$length <- gLength(rp, byid=TRUE) / 1000
x <- tapply(rp$length, rp$layer, sum)
r <- raster(rs)
r[as.integer(names(x))] <- x
Robert Hijmans
sumber
Tampaknya metode yang paling efisien dan elegan dari yang terdaftar. Juga, belum pernah melihat raster::intersect() sebelumnya, saya suka itu menggabungkan atribut fitur berpotongan, tidak seperti rgeos::gIntersection().
Matt SM
+1 selalu menyenangkan melihat solusi yang lebih efisien!
Jeffrey Evans
@ Robertt, saya mencoba menggunakan solusi ini untuk masalah lain, di mana saya ingin melakukan hal yang sama seperti yang ditanyakan di utas ini, tetapi dengan shapefile jalan yang sangat besar (untuk seluruh benua). Tampaknya telah bekerja, tetapi ketika saya mencoba melakukan gambar seperti yang dilakukan oleh @ fdetsch, saya tidak mendapatkan kotak yang berdekatan tetapi hanya beberapa kotak berwarna pada gambar (lihat di tinypic.com/r/20hu87k/9 ).
Doon_Bogan
Dan dengan paling efisien ... dengan dataset sampel saya: run time 0.6sec untuk solusi ini vs 8.25sec untuk solusi paling upvotes
user3386170
1

Anda tidak perlu for for loop. Hanya memotong semuanya sekaligus dan kemudian menambahkan panjang baris ke segmen baris baru menggunakan fungsi "SpatialLinesLengths" di sp. Kemudian, menggunakan fungsi raster paket paket raster dengan argumen fun = sum Anda bisa membuat raster dengan jumlah panjang garis yang memotong setiap sel. Menggunakan jawaban di atas dan data terkait di sini adalah kode yang akan menghasilkan hasil yang sama.

require(rgdal)
require(raster)
require(sp)
require(rgeos)

setwd("D:/TEST/RDSUM")
roads <- readOGR(getwd(), "TZA_roads")
  roads <- spTransform(roads, CRS("+init=epsg:21037"))
    rrst <- raster(extent(roads), crs=projection(roads))

# Intersect lines with raster "polygons" and add length to new lines segments
rrst.poly <- rasterToPolygons(rrst)
  rp <- gIntersection(roads, rrst.poly, byid=TRUE)
    rp <- SpatialLinesDataFrame(rp, data.frame(row.names=sapply(slot(rp, "lines"), 
                               function(x) slot(x, "ID")), ID=1:length(rp), 
                               length=SpatialLinesLengths(rp)/1000) ) 

# Rasterize using sum of intersected lines                            
rd.rst <- rasterize(rp, rrst, field="length", fun="sum")

# Plot results
require(RColorBrewer)
spplot(rd.rst, scales = list(draw=TRUE), xlab="x", ylab="y", 
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout=list("sp.lines", rp), 
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))
Jeffrey Evans
sumber
Pertama kali saya melihat SpatialLinesLengths. Kira tidak ada kata terlambat untuk belajar, terima kasih (: rasterizebutuh waktu cukup lama, (7 kali lebih lama dari pendekatan teratas di mesin saya)
fdetsch
Saya memang memperhatikan bahwa rasterisasi lambat. Namun, untuk masalah besar, loop for akan memperlambat segalanya dan saya pikir Anda akan melihat solusi yang jauh lebih optimal dengan rasterize. Selain itu pengembang paket raster berusaha keras untuk membuat setiap rilis lebih optimal dan lebih cepat.
Jeffrey Evans
Salah satu masalah potensial yang saya temukan dengan teknik ini adalah bahwa rasterize()fungsinya mencakup semua garis yang menyentuh sel tertentu. Ini menghasilkan panjang segmen garis yang dihitung dua kali dalam beberapa kasus: sekali di sel yang seharusnya dan di sel yang berdekatan yang disentuh titik akhir garis.
Matt SM
Ya, tetapi "rp" adalah objek yang dirasterisasi yang merupakan perpotongan poligon dan titik.
Jeffrey Evans
1

Inilah pendekatan lain. Itu menyimpang dari yang sudah diberikan dengan menggunakan spatstatpaket. Sejauh yang saya tahu, paket ini memiliki versi objek spasial sendiri (misalnya imvs rasterobjek), tetapi maptoolspaket ini memungkinkan konversi bolak-balik antara spatstatobjek dan objek spasial standar.

Pendekatan ini diambil dari pos R-sig-Geo ini .

require(sp)
require(raster)
require(rgdal)
require(spatstat)
require(maptools)
require(RColorBrewer)

# Load data and transform to UTM
roads <- shapefile('data/TZA_roads.shp')
roadsUTM <- spTransform(roads, CRS("+init=epsg:21037"))

# Need to convert to a line segment pattern object with maptools
roadsPSP <- as.psp(as(roadsUTM, 'SpatialLines'))

# Calculate lengths per cell
roadLengthIM <- pixellate.psp(roadsUTM, dimyx=10)

# Convert pixel image to raster in km
roadLength <- raster(dtanz / 1000, crs=projection(roadsUTM))

# Plot
spplot(rtanz, scales = list(draw=TRUE), xlab="x", ylab="y", 
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout=list("sp.lines", roadsUTM), 
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))

Imgur

Bit paling lambat adalah mengubah jalan dari SpatialLineske Pola Segmen Garis (yaitu spatstat::psp). Setelah selesai, bagian penghitungan panjang sebenarnya cukup cepat, bahkan untuk resolusi yang jauh lebih tinggi. Misalnya, di MacBook 2009 lama saya:

system.time(pixellate(tanzpsp, dimyx=10))
#    user  system elapsed 
#   0.007   0.001   0.007

system.time(pixellate(tanzpsp, dimyx=1000))
#    user  system elapsed 
#   0.146   0.032   0.178 
Matt SM
sumber
Hmmmm ... Saya sungguh berharap kapak itu tidak ada dalam notasi ilmiah. Adakah yang tahu cara memperbaikinya?
Matt SM
Anda dapat memodifikasi pengaturan R global dan mematikan notasi menggunakan: options (scipen = 999)) namun, saya tidak tahu apakah kisi akan menghormati pengaturan lingkungan global.
Jeffrey Evans
0

Biarkan saya hadir Anda paket vena dengan beberapa fungsi untuk bekerja garis spasial dan impor sf dan data.table

library(vein)
library(sf)
library(cptcity)
data(net)
netsf <- st_as_sf(net) #Convert Spatial to sf
netsf <- st_transform(netsf, 31983) # Project data
netsf$length_m  <- st_length(netsf)
netsf <- netsf[, "length_m"]
g <- make_grid(netsf, width = 1000) #Creat grid of 1000m spacing with columns id for each feature
# Number of lon points: 12
# Number of lat points: 11

gnet <- emis_grid(netsf, g)
plot(gnet["length_m"])

sf_to_raster <- function(x, column, ncol, nrow){
  x <- sf::as_Spatial(x)
  r <- raster::raster(ncol = ncol, nrow = nrow)
  raster::extent(r) <- raster::extent(x)
  r <- raster::rasterize(x, r, column)
  return(r)
}

rr <- sf_to_raster(gnet, "length_m", 12, 11)
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = cpt(), scales = list(draw = T))
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = cpt(pal = 5176), scales = list(draw = T))
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = lucky(), scales = list(draw = T))
# Colour gradient: neota_flor_apple_green, number: 6165

masukkan deskripsi gambar di sini

masukkan deskripsi gambar di sini

masukkan deskripsi gambar di sini

Sergio
sumber
-1

Ini mungkin terdengar agak naif tetapi jika itu adalah sistem jalan maka pilih jalan dan simpan ke papan klip kemudian lihat menemukan alat yang memungkinkan Anda untuk menambahkan buffer ke clipboard mengaturnya ke lebar jalan yang sah yaitu 3 meter +/- ingat bahwa penyangga adalah dari garis tengah ke tepi * 2 i untuk setiap sisi sehingga penyangga 3 meter sebenarnya adalah jalan 6 meter dari sisi ke sisi.

derekB
sumber
Apa hubungannya lebar jalan dengan panjang jalan? Jawaban ini tidak menjawab pertanyaan.
alphabetasoup