Szybsze przetwarzanie wektora na raster za pomocą R.

Przekształcam wektor na raster w R. Jednak proces ten był zbyt długi. Czy istnieje możliwość włączenia skryptu w proces wielowątkowy lub przetwarzanie GPU, aby zrobić to szybciej?

Mój skrypt do zrasteryzowanego wektora.

r.raster = raster()
extent(r.raster) = extent(setor) #definindo o extent do raster
res(r.raster) = 10 #definindo o tamanho do pixel
setor.r = rasterize(setor, r.raster, 'dens_imov')

r.raster

klasa: RasterLayer wymiary: 9636, 11476, 110582736 (nrow, ncol, ncell) rozdzielczość: 10, 10 (x, y) zasięg: 505755, 620515, 8555432, 8651792 (xmin, xmax, ymin, ymax) koordyn. ref. : + proj = longlat + datum = WGS84 + ellps = WGS84 + towgs84 = 0,0,0

setor

cechy klasy: SpatialPolygonsDataFrame: 5419 zasięg: 505755, 620515.4, 8555429, 8651792 (xmin, xmax, ymin, ymax) koordyn. ref. : + proj = utm + zone = 24 + południe + ellps = jednostki GRS80 + = m + zmienne no_defs: 6 nazw: ID, CD_GEOCODI, TIPO, dens_imov, area_m, domicilios 1 min wartości: 35464, 290110605000001, RURAL, 0,00000003,100004, 1,0000 wartości maksymalne: 58468, 293320820000042, URBANO, 0,54581673,99996, 99,0000

Wydruk setora

— Diogo Caribé
źródło

Czy możesz publikować streszczenia Setora i r.raster? Chciałbym mieć pojęcie o liczbie obiektów w zestawie i wymiarach r.raster. wystarczy wydrukować je w porządku

— mdsumner,

Podsumowanie zamieszczam w treści pytania.

— Diogo Caribé,

Nie podsumowanie, po prostu wydrukuj - informacje, o które prosiłem, a nie tgere

— mdsumner,

Przepraszam, umieszczam nadruk.

— Diogo Caribé,

Ach, rozczarowany, nie pomyślałem o tym, dopóki nie zobaczyłem wydruku - upewnij się, że rzut rastra pasuje do wielokątów, w tym momencie tak nie jest - spróbuj r <- raster (setor); res (r) <- 10; setor.r = rasterize (setor, r, 'dens_imov') - ale także spróbuj najpierw ustawić res (r) <- 250, żebyś miał pojęcie, jak długo potrwa wersja wysokiej rozdzielczości

— mdsumner

Próbowałem „zrównoleglić” funkcję rasterizeza pomocą Rpakietu parallelw następujący sposób:

podzielić obiekt SpatialPolygonsDataFrame na nczęści
rasterize każda część osobno
scal wszystkie części w jeden raster

W moim komputerze rasterizefunkcja zrównoleglona zajęła 2,75 razy mniej niż funkcja bez zrównoleglenia rasterize.

Uwaga: poniższy kod pobiera plik kształtu wielokąta (~ 26,2 MB) z Internetu. Możesz użyć dowolnego obiektu SpatialPolygonDataFrame. To tylko przykład.

Załaduj biblioteki i przykładowe dane:

# Load libraries
library('raster')
library('rgdal')

# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)

# Convert NAMES level 2 to factor 
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)

# Plot BRA_adm2
plot(BRA_adm2)
box()

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(BRA_adm2)

# Define pixel size
res(r.raster) <- 0.1

Rycina 1: Brazylia SpatialPolygonsDataFrame

Prosty przykład wątku

# Simple thread -----------------------------------------------------------

# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))

Czas na moim laptopie:

# Output:
# user  system elapsed 
# 23.883    0.010   23.891

Przykład wątku wielowątkowego

# Multithread -------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library('parallel')

# Calculate the number of cores
no_cores <- detectCores() - 1

# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])

# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)

# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))

# Finish
stopCluster(cl)

# Merge all raster parts
rMerge <- do.call(merge, rParts)

# Plot raster
plot(rMerge)

Rycina 2: Wykres rastrowy Brazylii

Czas na moim laptopie:

# Output:
# user  system elapsed 
# 0.203   0.033   8.688

Więcej informacji o równoległości w R :

— Guzmán
źródło

Bardzo dobra odpowiedź!

— Diogo Caribé,

Czy nie ustawiasz n jako liczby rdzeni na komputerze?

— Sam

@Sam Myślę, że powinno działać bez problemu, ale nie wiem, czy jest lepiej, czy nie! Założyłem, że jeśli podzielę funkcje na części n równe liczbie rdzeni, być może jedna z tych części mogłaby być łatwiejsza w obróbce, a rdzeń, który ją przetwarzał, byłby bezużyteczny! Jeśli jednak masz więcej części niż rdzeni, gdy jeden rdzeń zakończy przetwarzanie jednej części, wziąłby on inną część. Ale z pewnością nie jestem pewien! Jakakolwiek pomoc w tej sprawie będzie mile widziana.

— Guzmán

Dzisiaj wieczorem przeprowadzę kilka testów. Na małym pliku kształtu (około 25 km na 25 km), zrasteryzowanym do 50 m, nieznacznie poprawiono użycie n = 2,4 lub 8 w stosunku do n = 20, 30 lub nawet do 50. Dziś wieczorem poddam się bardzo dużemu plikowi kształtu i zrasteryzować do 25m. Przetwarzanie z jednym rdzeniem trwa 10 godzin, więc zobaczymy, co robią różne wartości n !! (n = 50 to niespełna 1 godzina)

— Sam

@ Guzmán Ponownie uruchamiam kod. Wystąpił jednak błąd i nie wiem dlaczego. Możesz mi pomóc? Błąd w checkForRemoteErrors (val): 7 węzłów wygenerowało błędy; pierwszy błąd: nie znaleziono obiektu „BRA_adm2”

— Diogo Caribé