Konwertuj plik kształtu linii na raster, wartość = całkowita długość linii w komórce

Mam plik kształtu linii reprezentujący sieć dróg. Chcę zrasteryzować te dane, a wartości wynikowe w rastrze pokazują całkowitą długość linii, które mieszczą się w komórce rastra.

Dane są w projekcji British National Grid, więc jednostkami będą metry.

Idealnie chciałbym wykonać tę operację przy użyciu Ri domyślam się, że rasterizefunkcja z rasterpakietu odegrałaby rolę w osiągnięciu tego, po prostu nie mogę ustalić, jaka powinna być zastosowana funkcja.

Po ostatnim pytaniu możesz skorzystać z funkcji oferowanych przez pakiet rgeos , aby rozwiązać problem. Ze względu na odtwarzalność pobrałem plik kształtu dróg tanzańskich z DIVA-GIS i umieściłem go w bieżącym katalogu roboczym. Do nadchodzących zadań potrzebne będą trzy pakiety:

• rgdal do ogólnego przetwarzania danych przestrzennych
• raster do rasteryzacji danych pliku shapefile
• rgeos, aby sprawdzić skrzyżowanie dróg z szablonem rastrowym i obliczyć długości dróg

W związku z tym Twoje pierwsze wiersze mogłyby wyglądać następująco:

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

Następnie musisz zaimportować dane pliku kształtu. Zauważ, że pliki kształtów DIVA-GIS są dystrybuowane w EPSG: 4326, więc wyświetlę plik kształtu do EPSG: 21037 (UTM 37S), aby zajmować się licznikami, a nie stopniami.

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

Do późniejszej rasteryzacji potrzebny będzie szablon rastrowy obejmujący zasięg przestrzenny pliku kształtu. Szablon rastrowy składa się domyślnie z 10 wierszy i 10 kolumn, co pozwala uniknąć zbyt długich czasów obliczeń.

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

Teraz, gdy szablon jest skonfigurowany, przejdź przez wszystkie komórki rastra (który obecnie składa się tylko z wartości NA). Przypisując wartość „1” do bieżącej komórki, a następnie wykonując rasterToPolygons, wynikowy plik kształtu „tmp_shp” automatycznie przechowuje zakres aktualnie przetwarzanego piksela. gIntersectswykrywa, czy ten zakres pokrywa się z drogami. Jeśli nie, funkcja zwróci wartość „0”. W przeciwnym razie plik kształtu drogi jest przycinany przez bieżącą komórkę, a całkowita długość „linii przestrzennych” w tej komórce jest obliczana za pomocą 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)
  }
})

Na koniec możesz wstawić obliczone długości (które są konwertowane na kilometry) do szablonu rastra i wizualnie zweryfikować swoje wyniki.

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

Poniżej zmodyfikowano rozwiązanie Jeffreya Evansa. To rozwiązanie jest znacznie szybsze, ponieważ nie używa rasteryzacji

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

Nie potrzebujesz pętli for. Po prostu przecinaj wszystko naraz, a następnie dodaj długości linii do nowych segmentów linii za pomocą funkcji „SpatialLinesLengths” w sp. Następnie, używając funkcji rasteryzacji pakietu rastrowego z argumentem fun = sum, możesz utworzyć raster z sumą długości linii przecinających każdą komórkę. Wykorzystanie powyższej odpowiedzi i powiązanych danych tutaj jest kodem, który wygeneruje te same wyniki.

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

Oto jeszcze inne podejście. Różni się od tych już podanych przy użyciu spatstatpakietu. O ile mi wiadomo, ten pakiet ma własną wersję obiektów przestrzennych (np. imVs. rasterobiekty), ale maptoolspakiet umożliwia konwersję tam iz powrotem między spatstatobiektami i standardowymi obiektami przestrzennymi.

To podejście zostało zaczerpnięte z tego postu 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))

Najwolniejszym bitem jest przekształcenie dróg z SpatialLineswzoru odcinka linii (tj spatstat::psp.). Po wykonaniu tej czynności obliczanie rzeczywistej długości jest dość szybkie, nawet w przypadku znacznie wyższych rozdzielczości. Na przykład na moim starym MacBooku 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 

Pozwól, że przedstawię Ci pakiet vein z kilkoma funkcjami do pracy z liniami przestrzennymi i importowania sf i 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

Może to zabrzmieć nieco naiwnie, ale jeśli jest to układ drogowy, wybierz drogi i zapisz w schowku, a następnie znajdź narzędzie, które pozwala dodać bufor do schowka, ustaw go na dozwoloną szerokość drogi, tj. 3 metry +/- pamiętaj, że bufor znajduje się od linii środkowej do krawędzi * 2 i dla każdej strony, więc 3-metrowy bufor jest w rzeczywistości 6-metrową drogą z boku na bok.