Convertir le fichier de formes de lignes en raster, valeur = longueur totale des lignes dans la cellule

J'ai un shapefile de ligne représentant un réseau routier. Je souhaite pixelliser ces données, les valeurs résultantes dans le raster indiquant la longueur totale des lignes qui se trouvent dans la cellule raster.

Les données sont dans la projection du British National Grid, les unités seront donc en mètres.

Idéalement, je voudrais effectuer cette opération en utilisant R, et je suppose que la rasterizefonction du rasterpackage jouerait un rôle dans la réalisation de cela, je ne peux tout simplement pas déterminer quelle devrait être la fonction appliquée.

Suite à une question récente , vous souhaiterez peut-être utiliser les fonctionnalités offertes par le package rgeos pour résoudre votre problème. Pour des raisons de reproductibilité, j'ai téléchargé un fichier de formes de routes tanzaniennes à partir de DIVA-GIS et l'ai mis dans mon répertoire de travail actuel. Pour les tâches à venir, vous aurez besoin de trois packages:

• rgdal pour la gestion générale des données spatiales
• raster pour la rastérisation des données du fichier de formes
• rgeos pour vérifier l'intersection des routes avec un modèle raster et calculer les longueurs des routes

Par conséquent, vos premières lignes de pourraient ressembler à ceci:

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

Après cela, vous devez importer les données du fichier de formes. Notez que les fichiers de formes DIVA-GIS sont distribués en EPSG: 4326, donc je projetterai le fichier de formes en EPSG: 21037 (UTM 37S) pour traiter les mètres plutôt que les degrés.

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

Pour la pixellisation ultérieure, vous aurez besoin d'un modèle de tramage qui couvre l'étendue spatiale de votre fichier de formes. Le modèle raster se compose de 10 lignes et 10 colonnes par défaut, évitant ainsi des temps de calcul trop longs.

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

Maintenant que le modèle est configuré, parcourez toutes les cellules du raster (qui se compose actuellement de valeurs NA uniquement). En attribuant une valeur de «1» à la cellule actuelle et en l'exécutant par la suite rasterToPolygons, le fichier de formes résultant «tmp_shp» contient automatiquement l'étendue du pixel actuellement traité. gIntersectsdétecte si cette étendue chevauche des routes. Sinon, la fonction renverra une valeur de «0». Sinon, le fichier de formes de route est rogné par la cellule actuelle et la longueur totale des «SpatialLines» dans cette cellule est calculée à l'aide de 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)
  }
})

Enfin, vous pouvez insérer les longueurs calculées (qui sont converties en kilomètres) dans le modèle raster et vérifier visuellement vos résultats.

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))

Ce qui suit est modifié à partir de la solution de Jeffrey Evans. Cette solution est beaucoup plus rapide car elle n'utilise pas la pixellisation

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

Vous n'avez pas besoin d'une boucle for. Il suffit de couper tout à la fois, puis d'ajouter des longueurs de ligne aux nouveaux segments de ligne à l'aide de la fonction "SpatialLinesLengths" de sp. Ensuite, en utilisant la fonction rasterize du package raster avec l'argument fun = sum, vous pouvez créer un raster avec la somme des longueurs de ligne coupant chaque cellule. L'utilisation de la réponse ci-dessus et des données associées ici est un code qui générera les mêmes résultats.

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))

Voici encore une autre approche. Il s'écarte de ceux déjà donnés en utilisant le spatstatpackage. Pour autant que je sache, ce package a sa propre version des objets spatiaux (par exemple, imvs rasterobjets), mais le maptoolspackage permet la conversion entre les spatstatobjets et les objets spatiaux standard.

Cette approche est tirée de ce poste 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))

Le bit le plus lent convertit les routes d' SpatialLinesun modèle de segment de ligne (c. spatstat::psp-à-d.). Une fois cela fait, les pièces de calcul de la longueur réelle sont assez rapides, même pour des résolutions beaucoup plus élevées. Par exemple, sur mon ancien 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 

Permettez-moi de vous présenter la veine du package avec plusieurs fonctions pour travailler les lignes spatiales et importer sf et 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

Cela peut sembler un peu naïf, mais s'il s'agit d'un système routier, sélectionnez les routes et enregistrez-les dans un presse-papiers, puis recherchez un outil qui vous permet d'ajouter un tampon au presse-papiers, réglez-le sur la largeur légale de la route, c'est-à-dire 3 mètres +/- rappelez-vous que le tampon est de la ligne médiane au bord * 2 i pour chaque côté, donc un tampon de 3 mètres est en fait une route de 6 mètres d'un côté à l'autre.