JeremyGelb / spNetwork

An R package to perform spatial analysis on networks.
GNU General Public License v2.0
35 stars 2 forks source link

Error in coordsLines[[nearest_line_index[x]]] #11

Closed sh3y1h closed 2 years ago

sh3y1h commented 2 years ago

Hi Jeremy,

Great work on spNetwork and thank you! I was trying to see how I could use the spNetwork to calculate oil spill density along the pipeline network, however, I keep getting stuck on this message:

Error in coordsLines[[nearest_line_index[x]]] : 
  attempt to select less than one element in get1index <real>

Here is my code and error message:

setwd("C:/Users/seyib/Downloads/OneDrive_2022-09-08/Getting confident with R/Dependencies")

pipeline_v2 <- readOGR(dsn = "C:/Users/seyib/Downloads/OneDrive_2022-09-08/Getting confident with R/Dependencies/Pipeline_test_run_qgis_reproject.shp")

spill_v3 <- readOGR(dsn = "C:/Users/seyib/Downloads/OneDrive_2022-09-08/Getting confident with R/Dependencies/spill_test_run_qgis_reproject.shp")

pipeline_v2_sf <- st_as_sf(pipeline_v2, coords = c("x", "y"), crs = 32632, agr = "constant")
st_agr(pipeline_v2_sf)

spill_v3_sf <- st_as_sf(spill_v3, coords = c("Longitude", "Latitude"), crs = 32632, agr = "constant")
st_agr(spill_v3_sf)

lixels <- lixelize_lines(pipeline_v2_sf,150, 50)
sample_points <- lines_center(lixels)
nkde_densities <- nkde(lines = pipeline_v2_sf,
                       events = spill_v3_sf,
                       w = rep(1,nrow(spill_v3_sf)),
                       samples = sample_points,
                       kernel_name = "quartic",
                       bw = 450,
                       adaptive = TRUE, trim_bw = 900,
                       method = "discontinuous",
                       div = "bw",
                       max_depth = 10,
                       digits = 2, tol = 0.1, agg = 5,
                       grid_shape = c(1,1),
                       verbose = FALSE)

spill_v3_sf (Event) pipeline_v2_sf (Network)

JeremyGelb commented 2 years ago

Hello There !

Thank you for your interest in spNetwork and the bug you reported. Could you send me the data your are using or a sample (also producing the bug) so that I could investigate the origin of the bug?

All the best!

sh3y1h commented 2 years ago

Thanks for your reply, Jeremy.

Please find attached the files (Network file: Pipeline (shp and extracted csv files); Events: spill (shp and extracted csv files).

Best regards, Seyi [Pipeline_test_run_qgis_reproject.zip](https://github.com/JeremyGelb/spNetwork/ [spill_test_run_qgis_reproject.zip](https://github.com/JeremyGelb/spNetwork/files/96027 pipeline_v2_sf.csv spill_v3_sf2.csv

23/spill_test_run_qgis_reproject.zip) files/9602721/Pipeline_test_run_qgis_reproject.zip)

JeremyGelb commented 2 years ago

Hello again,

I checked your case and the error does not come from the package but your data. If you take a quick look at them on Qgis you could notice this :

image

Some points are very far away from the network. This is causing the error.

I adapted your code to filter points that are too far (500 meters away from their closest line). Also, your study area is huge, you should consider using larger lixels (150m is too small for a study area larger than 600km), otherwise, mapping will be difficult. Finally, it seems that some parts of your network are not truly connected, here is an example: image

Please, check if this is normal or an error in the data. In that case, you should edit the features in a GIS software (like Qgis).

Here is an example of an adapted code, but keeping your parameters and working with the current version of spNetwork on github.

library(spNetwork)
library(sf)
library(tmap)

# Loading the data with the package sf (no need of readOGR)
setwd("C:/Users/Gelb/Desktop/TEMP/spNetwork_issue12")
pipeline_v2 <- st_read("Pipeline_test_run_qgis_reproject.shp")
spill_v3 <- st_read("spill_test_run_qgis_reproject.shp")

# Creating the lixels along the lines (NOTE : 150m is probably too small)
lixels <- lixelize_lines(pipeline_v2,150, 50)

sample_points <- lines_center(lixels)

# Filtering the points that are 500 meters away from a line
distances <- st_distance(spill_v3, pipeline_v2)
Best_dist <- apply(distances, 1, min)
ok_pts <- subset(spill_v3, Best_dist < 500)

# Calculating densities
nkde_densities <- nkde(lines = pipeline_v2,
                       events = ok_pts,
                       w = rep(1,nrow(ok_pts)),
                       samples = sample_points,
                       kernel_name = "quartic",
                       bw = 450,
                       adaptive = TRUE, trim_bw = 900,
                       method = "discontinuous",
                       div = "bw",
                       max_depth = 10,
                       digits = 2, tol = 0.1, agg = 5,
                       grid_shape = c(1,1),
                       verbose = TRUE)

# quick mapping (slow because of the lixels' size)
lixels$density <- nkde_densities$k
tm_shape(lixels) +
  tm_lines("density", n = 7, style = "jenks", palette = "viridis")
sh3y1h commented 2 years ago

Thanks, Jeremy 👍 !

Except for the final line, everything appears to be in order now:

tm_lines(col = "density", n = 7, style = "jenks", palette = "viridis")

It's returning this error:

_Error: Invalid line colors

However, if I change "Density" to "NA," it runs, but the printed map is plain lines with no density. For now, I'm moving on to the Spatio-temporal part of the package, as my data is temporal.

sh3y1h commented 2 years ago

Hello again, Jeremy,

The time you've spent attending to this has been much appreciated. Attempting the Spatio-temporal section, I couldn't go forward because of some errors; below are the lines of code and the resulting errors encountered while calculating Spatio-temporal.

# To remove rows with NAs
ok_pts <- na.omit(ok_pts)

convert the time to numeric

ok_pts$Time <- as.POSIXct(ok_pts$Incident_D, format = "%Y/%m/%d") start <- as.POSIXct("2019/01/02", format = "%Y/%m/%d")

ok_pts$Time <- difftime(ok_pts$Incident_D, start, units = "days") ok_pts$Time <- as.numeric(ok_pts$Incident_D)

months <- as.character(1:12) months <- ifelse(nchar(months)==1, paste0("0", months), months) months_starts_labs <- paste("2019/",months,"/05", sep = "") months_starts_num <- as.POSIXct(months_starts_labs, format = "%Y/%m/%d") months_starts_num <- difftime(months_starts_num, start, units = "days") months_starts_num <- as.numeric(months_starts_num) months_starts_labs <- gsub("2019/", "", months_starts_labs, fixed = TRUE)

``

choosing sample in times (every 10 days)

sample_time <- seq(0, max(ok_pts$Time), 10)

Error: Error in seq.default(0, max(ok_pts$Time), 10) : 'to' must be a finite number In addition: Warning message: In max(ok_pts$Time) : no non-missing arguments to max; returning -Inf

calculating densities

tnkde_densities <- tnkde(lines = pipeline_v2, events = ok_pts, time_field = "Incident_D", w = rep(1, nrow(ok_pts)), samples_loc = sample_points, samples_time = sample_time, kernel_name = "quartic", bw_net = 700, bw_time = 60, adaptive = TRUE, trim_bw_net = 900, trim_bw_time = 80, method = "discontinuous", div = "bw", max_depth = 10, digits = 2, tol = 0.01, agg = 15, grid_shape = c(1,1), verbose = FALSE)


Error: Error in s - events : non-numeric argument to binary operator
JeremyGelb commented 2 years ago

Hello again !

Sorry but this is not exactly a place to be helped with general R questions. I will give you some help but I do not have the time to investigate errors caused by functions that are not related to spNetwork. Don't get me wrong, it is great that you are exploring R and working with this amazing software! I am also very happy that you are using spNetwork in your work. But the (my free) time I am spending on this type of issue cannot be spent on improving the package.

So let us see what is wrong here.

ok_pts$Time <- as.POSIXct(ok_pts$Incident_D, format = "%Y/%m/%d")

This cannot work because of your date format : "30/05/2019" (in the dataframe ok_pts) Look : "%Y/%m/%d" means YEAR/MONTH/DAY, and your data are "DAY/MONTH/YEAR" So change the parameter format for : "%d/%m/%Y"

ok_pts$Time <- difftime(ok_pts$Incident_D, start, units = "days")

Here, you are using the column Incident_D in the function difftime. Incident_D is a character column, you must use the DateTime column you just created (Time), with the function as.POSIXct.

ok_pts$Time_num <- difftime(ok_pts$Time, start, units = "days")
ok_pts$Time_num <- as.numeric(ok_pts$Time_num)
sample_time <- seq(0, max(ok_pts$Time_num), 10)

And then, I will works like a charm!

Here is the full corrected code.

library(spNetwork)
library(sf)
library(tmap)

setwd("C:/Users/Gelb/Desktop/TEMP/spNetwork_issue12")

pipeline_v2 <- st_read("Pipeline_test_run_qgis_reproject.shp")

spill_v3 <- st_read("spill_test_run_qgis_reproject.shp")

lixels <- lixelize_lines(pipeline_v2,150, 50)

sample_points <- lines_center(lixels)

distances <- st_distance(spill_v3, pipeline_v2)
Best_dist <- apply(distances, 1, min)

ok_pts <- subset(spill_v3, Best_dist < 500)

# To remove rows with NAs
ok_pts <- na.omit(ok_pts)

# creating a numeric vector for temporal column
ok_pts$Time <- as.POSIXct(ok_pts$Incident_D, format = "%d/%m/%Y")
start <- as.POSIXct("2019/01/02", format = "%Y/%m/%d")

ok_pts$Time_num <- difftime(ok_pts$Time, start, units = "days")
ok_pts$Time_num <- as.numeric(ok_pts$Time_num)

sample_time <- seq(0, max(ok_pts$Time_num), 10)

tnkde_densities <- tnkde(lines = pipeline_v2,
                       events = ok_pts,
                       time_field = "Time_num",
                       samples_loc = sample_points,
                       samples_time = sample_time,
                       w = rep(1,nrow(ok_pts)),
                       kernel_name = "quartic",
                       bw_net = 700, bw_time = 60,
                       adaptive = TRUE,
                       trim_bw_net = 900,
                       trim_bw_time = 80,
                       method = "discontinuous",
                       div = "bw",
                       max_depth = 10,
                       digits = 2, tol = 0.1, agg = 5,
                       grid_shape = c(1,1),
                       verbose = TRUE)

Maybe the issue here is the lack of documentation? Have you worked with the vignettes and the examples in the package? Have you checked the website for more examples?

All the best !

sh3y1h commented 2 years ago

Thank you so much, Jeremy!

First and foremost, please accept my apologies for interrupting your valuable time. As you may have observed, I'm still new to R; thank you very much for taking the time to examine the code and for your patience. It worked like magic, as you stated! Finally, many thanks for the fantastic spNetwork package!

JeremyGelb commented 2 years ago

No problem, I am happy that I could help you! Keep exploring and have fun with R!

sh3y1h commented 2 years ago

Thanks for your reply, Jeremy.

Please find attached the files (Network file: Pipeline (shp and extracted csv files); Events: spill (shp and extracted csv files).

Best regards, Seyi

On Mon, 19 Sept 2022 at 19:42, JeremyGelb @.***> wrote:

Hello There !

Thank you for your interest in spNetwork and the bug you reported. Could you send me the data your are using or a sample (also producing the bug) so that I could investigate the origin of the bug?

All the best!

— Reply to this email directly, view it on GitHub https://github.com/JeremyGelb/spNetwork/issues/11#issuecomment-1251405842, or unsubscribe https://github.com/notifications/unsubscribe-auth/A3C2L5ZHWAAAXQ2IJ72OLF3V7CXYPANCNFSM6AAAAAAQQJWK2I . You are receiving this because you authored the thread.Message ID: @.***>