Оптимизация кода для подсчета очков внутри буферного кольца - PullRequest
Новая
       12

Оптимизация кода для подсчета очков внутри буферного кольца

0 голосов
/ 26 февраля 2020

Я использую модель эпидемии c, и у меня есть код ниже, чтобы сделать это. В целом, код выполняет следующие действия: 1 - создает растр с заданным разрешением, охватывающим расширение территории Бразилии; 2. Создайте точки сканирования из центроидов растровых пикселей и исключите точки из Бразилии; 3-Создать буферное кольцо из каждой точки сканирования с заданным радиусом; 4-Подсчитайте, сколько точек появилось за каждый год и каждую точку сканирования, и сохраните ее в матрице результатов.

Моя проблема в том, что я должен выполнить одну и ту же процедуру для многих разрешений растра и множества радиусов буферных колец. Когда я уменьшу разрешение, будет больше растровых точек для анализа, поэтому для выполнения кода требуется слишком много времени. Есть ли лучший способ сделать то же самое, но с гораздо более высокой производительностью?

Точки данных для подсчета находятся по этой ссылке: Точки данных 

# Cleaning the work environment
rm(list = ls())

# Calling packages
library(rgdal)
library(raster)
require(sf)
require(maptools)
require(geosphere)
library(rgeos)

# Getting shapefile from Brazil
bra <- getData("GADM", country="BRA", level=1)

# Setting a coordination reference system-CRS
bcrs <- "+proj=poly +lat_0=0 +lon_0=-54 +x_0=5000000 +y_0=10000000 +ellps=GRS80"

# Setting the CRS to Brazil
pbra <- sp::spTransform(bra, bcrs)

# Reading prosumers data
pros <- as(st_read("~/points.shp"), 'Spatial')
pros <- sp::spTransform(pros, bcrs)

# Creating a raster res=the number in meters of the cell size
r <- raster(pbra, res=1000000)

# Coordinates of raster centroids
scan <- SpatialPoints(coordinates(r), proj4string=CRS(bcrs))

# Extract points out of Brazil
scan_subset <- scan[pbra, ]

# Create buffer rings from scan points
ring <- st_as_sf(gBuffer(scan_subset, width=100000, byid=TRUE))

# Tranform to sf object
pros <- st_as_sf(pros)
pros$Data_Conex <- as.Date(pros$Data_Conex, "%d/%m/%Y")

# Create a result marix
result <- matrix(ncol = 3, nrow = dim(ring)[1]*8)
k=1
for (i in 1:dim(ring)[1]) {
    result[k:(k+7),1] <- i
    result[k:(k+7),2] <- 2012:2019
    k=k+8
}
# Start the clock!
ptm <- proc.time()
for (i in 1:dim(result)[1]) {
    result[i,3] <- dim(st_intersection(x = ring[result[i,1],1],
                                       y = subset(pros, 
                                                  pros$Data_Conex >= as.Date(paste0(result[i,2],"-01-01")) 
                                                  & pros$Data_Conex <= as.Date(paste0(result[i,2],"-12-31")))))[1]
}
# Stop the clock
proc.time() - ptm
...