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

14

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.

JPD
źródło
Może vignette('over', package = 'sp')może pomóc.

Odpowiedzi:

12

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

lines2raster

fdetsch
źródło
To jest niesamowite! Zmieniłem sapply()się pbsapply()i użył argumentu klastra cl = detectCores()-1. Teraz mogę uruchomić ten przykład równolegle!
philiporlando
11

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
Robert Hijmans
źródło
Wydaje się najbardziej wydajną i elegancką metodą spośród wymienionych. Ponadto, nie widziałem raster::intersect() wcześniej, podoba mi się, że łączy atrybuty przecinających się elementów, w przeciwieństwie do rgeos::gIntersection().
Matt SM
+1 zawsze miło widzi bardziej wydajne rozwiązania!
Jeffrey Evans,
@RobertH, próbowałem użyć tego rozwiązania dla innego problemu, w którym chcę zrobić to samo, co zadałem w tym wątku, ale z bardzo dużym plikiem kształtów dróg (dla całego kontynentu). Wygląda na to, że zadziałało, ale kiedy próbuję wykonać figurę wykonaną przez @ fdetsch, nie widzę ciągłej siatki, ale tylko kilka kolorowych kwadratów na zdjęciu (patrz w tinypic.com/r/20hu87k/9 ).
Doon_Bogan
I przez najbardziej wydajny ... z moim przykładowym zestawem danych: czas działania 0,6 s dla tego rozwiązania w porównaniu z 8,25 s dla najbardziej pozytywnego rozwiązania.
user3386170,
1

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))
Jeffrey Evans
źródło
Pierwszy raz widzę SpatialLinesLengths. Chyba nigdy nie jest za późno na naukę, dziękuję (: rasterizetrwa to jednak dość długo (7 razy dłużej niż górne podejście na moim komputerze).
fdetsch,
Zauważyłem, że rasteryzacja przebiega powoli. Jednak w przypadku dużych problemów pętla for naprawdę spowolni i myślę, że zobaczysz znacznie bardziej zoptymalizowane rozwiązanie dzięki rasterize. Poza tym twórca pakietu rastrowego stara się, aby każde wydanie było bardziej zoptymalizowane i szybsze.
Jeffrey Evans,
Jednym z potencjalnych problemów, które znalazłem przy tej technice, jest to, że rasterize()funkcja obejmuje wszystkie linie dotykające danej komórki. W niektórych przypadkach długości segmentów linii są liczone dwukrotnie: raz w komórce, którą powinni, i raz w sąsiedniej komórce, której dotyka punkt końcowy linii.
Matt SM
Tak, ale „rp” to obiekt rasteryzowany, który jest przecięciem wielokątów i punktów.
Jeffrey Evans
1

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

Imgur

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 
Matt SM
źródło
Hmmmm ... z pewnością chciałbym, żeby te osie nie były w notacji naukowej. Czy ktoś ma pomysł, jak to naprawić?
Matt SM
Możesz zmodyfikować globalne ustawienia R i wyłączyć notację używając: options (scipen = 999)), jednak nie wiem, czy sieć będzie honorować globalne ustawienia środowiska.
Jeffrey Evans
0

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

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

Sergio
źródło
-1

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.

derekB
źródło
Co szerokość drogi ma wspólnego z długością drogi. Ta odpowiedź nie dotyczy pytania.
alphabetasoup