Converter arquivo de forma de linha em raster, valor = comprimento total de linhas na célula

14

Eu tenho um shapefile de linha representando uma rede de estradas. Desejo rasterizar esses dados, com os valores resultantes na varredura mostrando o comprimento total das linhas que se enquadram na célula de varredura.

Os dados estão na projeção da British National Grid, portanto as unidades serão metros.

Idealmente, eu gostaria de executar esta operação usando R, e estou supondo que a rasterizefunção do rasterpacote tenha um papel importante nesse processo, simplesmente não consigo descobrir qual deve ser a função aplicada.

JPD
fonte
Talvez vignette('over', package = 'sp')possa ajudar.

Respostas:

12

Após uma pergunta recente , convém usar as funcionalidades oferecidas pelo pacote rgeos para resolver seu problema. Por motivos de reprodutibilidade, baixei um arquivo de forma das estradas da Tanzânia do DIVA-GIS e o coloquei no meu diretório de trabalho atual. Para as próximas tarefas, você precisará de três pacotes:

  • rgdal para manipulação geral de dados espaciais
  • rasterização para rasterização dos dados do shapefile
  • rgeos para verificar a interseção de estradas com modelo raster e calcular comprimentos de estrada

Conseqüentemente, suas primeiras linhas de could devem ficar assim:

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

Depois disso, você precisa importar os dados do shapefile. Observe que os shapefiles do DIVA-GIS são distribuídos no EPSG: 4326, então projetarei o shapefile no EPSG: 21037 (UTM 37S) para lidar com medidores e não com graus.

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

Para rasterização subsequente, você precisará de um modelo raster que cubra a extensão espacial do seu shapefile. O modelo de varredura consiste em 10 linhas e 10 colunas por padrão, evitando assim tempos de computação muito extensos.

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

Agora que o modelo está configurado, percorra todas as células da varredura (que atualmente consiste apenas em valores de NA). Ao atribuir um valor de '1' à célula atual e subsequentemente executar rasterToPolygons, o shapefile resultante 'tmp_shp' mantém automaticamente a extensão do pixel processado no momento. gIntersectsdetecta se essa extensão se sobrepõe às estradas. Caso contrário, a função retornará um valor '0'. Caso contrário, o shapefile da estrada é cortado pela célula atual e o comprimento total de 'SpatialLines' nessa célula está sendo calculado usando 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)
  }
})

Por fim, você pode inserir os comprimentos calculados (que são convertidos em quilômetros) no modelo de varredura e verificar visualmente seus resultados.

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
fonte
Isso é incrível! Mudei sapply()para pbsapply()e usei o argumento cluster cl = detectCores()-1. Agora eu sou capaz de executar este exemplo em paralelo!
Philiporlando
11

O abaixo é modificado a partir da solução de Jeffrey Evans. Esta solução é muito mais rápida, pois não utiliza rasterização

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
fonte
Parece o método mais eficiente e elegante dos listados. Além disso, nunca tinha visto raster::intersect() antes, gosto que ele combina os atributos dos recursos cruzados, ao contrário rgeos::gIntersection().
Matt SM
+1 sempre bom ver soluções mais eficientes!
Jeffrey Evans
@ RobertH, tentei usar esta solução para outro problema, onde quero fazer o mesmo que solicitado neste tópico, mas com um grande formato de arquivo de estradas (para todo um continente). Parece ter funcionado, mas quando tento fazer a figura como @ fdetsch, não recebo uma grade contígua, mas apenas alguns quadrados coloridos na imagem (veja em tinypic.com/r/20hu87k/9 ).
Doon_Bogan 23/05
E o mais eficiente ... com o meu conjunto de dados de amostra: tempo de execução 0,6seg para esta solução vs. 8,25s para a solução mais positiva.
user3386170
1

Você não precisa de um loop for. Basta interceptar tudo de uma só vez e adicionar comprimentos de linha aos novos segmentos de linha usando a função "SpatialLinesLengths" em sp. Em seguida, usando a função rasterize do pacote raster com o argumento fun = sum, é possível criar uma varredura com a soma do (s) comprimento (s) da linha que cruzam cada célula. Usando a resposta acima e os dados associados, aqui está o código que irá gerar os mesmos resultados.

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
fonte
Primeira vez que vejo SpatialLinesLengths. Acho que nunca é tarde demais para aprender, obrigado (: rasterizeleva muito longo, embora (7 vezes mais do que a abordagem superior em minha máquina).
fdetsch
Eu notei que a rasterização estava sendo lenta. No entanto, para problemas grandes, um loop for vai realmente desacelerar as coisas e acho que você verá uma solução muito mais otimizada com rasterize. Além disso, o desenvolvedor do pacote raster está se esforçando para tornar cada versão mais otimizada e mais rápida.
Jeffrey Evans
Um possível problema que encontrei com essa técnica é que a rasterize()função inclui todas as linhas que tocam uma determinada célula. Isso resulta em comprimentos de segmentos de linha sendo contados duas vezes em alguns casos: uma vez na célula eles deveriam e uma vez na célula adjacente que apenas o ponto final da linha toca.
Matt SM
Sim, mas "rp" é o objeto que está sendo rasterizado, que é a interseção de polígonos e pontos.
Jeffrey Evans
1

Aqui está mais uma abordagem. Ele se desvia daqueles já fornecidos usando o spatstatpacote. Tanto quanto posso dizer, este pacote tem sua própria versão de objetos espaciais (por exemplo, imvs. rasterobjetos), mas o maptoolspacote permite a conversão entre spatstatobjetos e objetos espaciais padrão.

Essa abordagem é retirada deste post do R-sig-Geo .

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

O bit mais lento é converter as estradas de SpatialLinespara um padrão de segmento de linha (ou seja spatstat::psp). Feito isso, as peças de cálculo do comprimento real são bastante rápidas, mesmo para resoluções muito mais altas. Por exemplo, no meu antigo MacBook de 2009:

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
fonte
Hmmmm ... eu certamente gostaria que esses eixos não estivessem em notação científica. Alguém tem alguma idéia de como consertar isso?
Matt SM
Você pode modificar as configurações globais de R e desativar a notação usando: options (scipen = 999)), no entanto, não sei se a estrutura atenderá às configurações globais do ambiente.
Jeffrey Evans
0

Deixe-me apresentar a veia do pacote com várias funções para trabalhar linhas espaciais e importar sf e 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

insira a descrição da imagem aqui

insira a descrição da imagem aqui

insira a descrição da imagem aqui

Sergio
fonte
-1

Isso pode parecer um pouco ingênuo, mas se for um sistema de estradas, selecione as estradas e salve em uma área de transferência. Procure uma ferramenta que permita adicionar um buffer à área de transferência, defina-a na largura legal da estrada, ou seja, 3 metros +/- lembre-se de que o buffer é da linha central até a borda * 2i para cada lado, de modo que um buffer de 3 metros é na verdade uma estrada de 6 metros de um lado para o outro.

derekB
fonte
O que a largura da estrada tem a ver com o comprimento da estrada. Esta resposta não aborda a questão.
alphabetasoup