Convierta el archivo de forma de línea en ráster, valor = longitud total de líneas dentro de la celda

14

Tengo un archivo de forma de línea que representa una red de carreteras. Deseo rasterizar estos datos, con los valores resultantes en el ráster que muestran la longitud total de las líneas que se encuentran dentro de la celda ráster.

Los datos se encuentran en la proyección de la National National Grid, por lo que las unidades serán metros.

Idealmente, me gustaría realizar esta operación utilizando R, y supongo que la rasterizefunción del rasterpaquete jugaría un papel en lograr esto, simplemente no puedo determinar cuál debería ser la función aplicada.

JPD
fuente
Quizás vignette('over', package = 'sp')pueda ayudar.

Respuestas:

12

Después de una pregunta reciente , es posible que desee utilizar las funcionalidades que ofrece el paquete rgeos para resolver su problema. Por razones de reproducibilidad, descargué un archivo shape de carreteras de Tanzania desde DIVA-GIS y lo puse en mi directorio de trabajo actual. Para las próximas tareas, necesitará tres paquetes:

  • rgdal para manejo de datos espaciales generales
  • ráster para rasterizar los datos del archivo de forma
  • rgeos para verificar la intersección de carreteras con plantilla ráster y calcular longitudes de carretera

En consecuencia, sus primeras líneas de podría debería verse así:

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

Después de eso, debe importar los datos del archivo de forma. Tenga en cuenta que los archivos de forma DIVA-GIS se distribuyen en EPSG: 4326, por lo que proyectaré el archivo de forma a EPSG: 21037 (UTM 37S) para tratar metros en lugar de grados.

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

Para la posterior rasterización, necesitará una plantilla de trama que cubra la extensión espacial de su shapefile. La plantilla ráster consta de 10 filas y 10 columnas de forma predeterminada, lo que evita tiempos de cálculo demasiado extensos.

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

Ahora que la plantilla está configurada, recorra todas las celdas del ráster (que actualmente consta solo de valores de NA). Al asignar un valor de '1' a la celda actual y luego ejecutarlo rasterToPolygons, el archivo de forma resultante 'tmp_shp' mantiene automáticamente la extensión del píxel procesado actualmente. gIntersectsdetecta si esta extensión se superpone con las carreteras. Si no, la función devolverá un valor de '0'. De lo contrario, la celda actual recorta el archivo de forma de la carretera y se calcula usando la longitud total de 'SpatialLines' dentro de esa celda 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)
  }
})

Finalmente, puede insertar las longitudes calculadas (que se convierten en kilómetros) en la plantilla ráster y verificar visualmente sus 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
fuente
¡Esto es asombroso! He cambiado sapply()a pbsapply()y se utiliza el argumento de clúster cl = detectCores()-1. ¡Ahora puedo ejecutar este ejemplo en paralelo!
philiporlando
11

Lo siguiente está modificado de la solución de Jeffrey Evans. Esta solución es mucho más rápida ya que no utiliza rasterizar

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
fuente
Parece el método más eficiente y elegante de los enumerados. Además, no lo había visto raster::intersect() antes, me gusta que combine los atributos de las características intersectadas, a diferencia rgeos::gIntersection().
Matt SM
¡+1 siempre es agradable ver soluciones más eficientes!
Jeffrey Evans
@RobertH, intenté usar esta solución para otro problema, donde quiero hacer lo mismo que se me pidió en este hilo, pero con un gran archivo de formas de carreteras (para todo un continente). Parece que funcionó, pero cuando trato de hacer la figura como lo hizo @ fdetsch, no obtengo una cuadrícula contigua sino solo unos pocos cuadrados de colores en la imagen (ver en tinypic.com/r/20hu87k/9 ).
Doon_Bogan
Y de la manera más eficiente ... con mi conjunto de datos de muestra: tiempo de ejecución de 0.6 segundos para esta solución frente a 8.25 segundos para la solución con más votos.
user3386170
1

No necesita un bucle for. Simplemente intersecte todo de una vez y luego agregue longitudes de línea a los nuevos segmentos de línea utilizando la función "SpatialLinesLengths" en sp. Luego, utilizando la función rasterizar paquete de ráster con el argumento fun = sum, puede crear un ráster con la suma de las longitudes de línea que intersecan cada celda. El uso de la respuesta anterior y los datos asociados aquí es un código que generará los mismos 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
fuente
Primera vez que veo SpatialLinesLengths. Supongo que nunca es demasiado tarde para aprender, gracias (: rasterizelleva bastante tiempo, aunque (7 veces más que el enfoque superior en mi máquina).
fdetsch
Me di cuenta de que rasterizar estaba siendo lento. Sin embargo, para problemas grandes, un bucle for realmente ralentizará las cosas y creo que verá una solución mucho más optimizada con rasterizar. Además, el desarrollador del paquete ráster se esfuerza por hacer que cada versión sea más optimizada y más rápida.
Jeffrey Evans
Un problema potencial que encontré con esta técnica es que la rasterize()función incluye todas las líneas que tocan una celda determinada. Esto da como resultado que las longitudes de segmento de línea se cuenten dos veces en algunos casos: una vez en la celda que se supone que deben hacer y una vez en la celda adyacente que el punto final de la línea solo toca.
Matt SM
Sí, pero "rp" es el objeto que se rasteriza, que es la intersección de polígonos y puntos.
Jeffrey Evans
1

Aquí hay otro enfoque. Se desvía de los que ya se proporcionan al usar el spatstatpaquete. Por lo que puedo decir, este paquete tiene su propia versión de objetos espaciales (por ejemplo, imcontra rasterobjetos), pero el maptoolspaquete permite la conversión de ida y vuelta entre spatstatobjetos y objetos espaciales estándar.

Este enfoque se toma de esta publicación 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

El bit más lento es convertir las carreteras de SpatialLinesa un patrón de segmento de línea (es decir spatstat::psp). Una vez hecho esto, las partes del cálculo de la longitud real son bastante rápidas, incluso para resoluciones mucho más altas. Por ejemplo, en mi vieja MacBook 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
fuente
Hmmmm ... desearía que esos ejes no estuvieran en notación científica. Alguien tiene alguna idea de cómo solucionar eso?
Matt SM
Puede modificar la configuración global de R y desactivar la notación usando: opciones (scipen = 999)) sin embargo, no sé si la red respetará la configuración del entorno global.
Jeffrey Evans
0

Permítame presentarle la veta del paquete con varias funciones para trabajar líneas espaciales e importar sf y 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

ingrese la descripción de la imagen aquí

ingrese la descripción de la imagen aquí

ingrese la descripción de la imagen aquí

Sergio
fuente
-1

Esto puede sonar un poco ingenuo, pero si se trata de un sistema de carreteras, seleccione las carreteras y guárdelas en un portapapeles, luego busque una herramienta que le permita agregar un búfer al portapapeles, ajústelo al ancho legal de la carretera, es decir, 3 metros +/- recuerde que el buffer es desde la línea central hasta el borde * 2 i para cada lado, por lo que un buffer de 3 metros es en realidad una carretera de 6 metros de lado a lado.

derekB
fuente
¿Qué tiene que ver el ancho del camino con la longitud del camino? Esta respuesta no aborda la pregunta.
alphabetasoup