0

I have a text file of coordinate data and I want to find the distance moved within a particular zone from that data. How would I go about that?

Here is a small sample of the data...

Time         X     Y     Z
============   ===== ===== =====
000:00:00.05    8.50 10.00  6.50
000:00:00.10    4.00 10.00  6.50
000:00:00.15    8.50 10.00  6.50
000:00:00.20    3.50 10.00  6.50
000:00:00.25    3.50 10.00  6.50
000:00:00.30    3.00 10.00  6.50
000:00:00.35    3.00 10.00  6.50
000:00:00.40    3.00 10.00  6.50
000:00:00.45    3.50 10.00  0.00
000:00:00.50    4.50 10.50  0.00
000:00:00.55    4.50 10.50  0.00
000:00:00.60    4.50 10.50  0.00
000:00:00.65    4.50 10.50  0.00
000:00:00.70    4.50 10.50  0.00
000:00:00.75    4.50 10.50  0.00
000:00:00.80    4.50 10.50  0.00
000:00:00.85    4.50 10.50  0.00
000:00:00.90    5.00 10.50  0.00
000:00:00.95    5.50 10.50  0.00
000:00:01.00    5.50 10.00  0.00
000:00:01.05    5.50 10.50  0.00
000:00:01.10    5.50 10.00  0.00
000:00:01.15    5.50 10.00  0.00
000:00:01.20    5.50 10.00  0.00
000:00:01.25    5.50 10.00  0.00
000:00:01.30    6.00 10.00  0.00
000:00:01.35    6.00 10.00  0.00
000:00:01.40    6.50 10.00  0.00
000:00:01.45    6.50  9.50  0.00
000:00:01.50    6.50  9.50  0.00
000:00:01.55    6.00  9.50  0.00
000:00:01.60    6.00  9.50  0.00
000:00:01.65    7.00  9.00  0.00
000:00:01.70    7.00  9.00  0.00
000:00:01.75    7.00  9.00  0.00
000:00:01.80    7.00  9.00  0.00
000:00:01.85    7.00  9.00  0.00
000:00:01.90    7.50  9.00  0.00
000:00:01.95    8.00  9.00  0.00
000:00:02.00    8.50  9.00  0.00
000:00:02.05    8.50  9.00  0.00
000:00:02.10    8.50  9.00  0.00
000:00:02.15    8.50  9.00  0.00
000:00:02.20    8.50  9.00  0.00
000:00:02.25    8.50  9.00  0.00
000:00:02.30    9.00  9.00  0.00
000:00:02.35    9.50  9.00  0.00
000:00:02.40    9.50  8.50  0.00
000:00:02.45    9.50  8.00  0.00
000:00:02.50   10.00  7.50  0.00
000:00:02.55   10.00  8.00  0.00
000:00:02.60   10.50  8.00  0.00
000:00:02.65   10.50  8.00  0.00
000:00:02.70   11.00  8.00  0.00
000:00:02.75   11.50  8.00  0.00
000:00:02.80   11.50  8.00  0.00
000:00:02.85   12.00  8.00  0.00
000:00:02.90   12.00  8.00  0.00
000:00:02.95   12.00  8.00  0.00
000:00:03.00   12.00  8.00  0.00
000:00:03.05   12.00  8.00  0.00
000:00:03.10   12.00  8.00  0.00
000:00:03.15   12.00  8.00  0.00
000:00:03.20   12.00  8.00  0.00
000:00:03.25   12.00  8.00  0.00
000:00:03.30   12.00  8.00  0.00
000:00:03.35   12.00  8.00  0.00
000:00:03.40   12.00  8.00  0.00
000:00:03.45   12.00  8.00  0.00
000:00:03.50   12.00  8.00  0.00
000:00:03.55   12.00  8.00  0.00
000:00:03.60   12.00  8.00  0.00
000:00:03.65   12.00  8.00  0.00
000:00:03.70   12.00  8.00  0.00
000:00:03.75   12.00  8.50  0.00
000:00:03.80   12.00  8.50  0.00
000:00:03.85   12.00  8.50  0.00
000:00:03.90   12.00  9.00  0.00
000:00:03.95   12.00  9.00  0.00
000:00:04.00   12.00  9.00  0.00

Basically, I have a 17 X 17 inch box and would like to measure the movement within the box. I would like to measure the center zone (maybe a 4 x 4 inch square) as well as the residual area zone. Each X and Y coordinate marks the 1 inch or 1//2 mark. You can ignore any Z coordinates.

Also, If I wanted to bin the movement by a selection of time (like 1 minute or so many seconds) how would I do that?

Here I fixed the time column to show it in second.milliseconds...

time     x     y     z
1    0.05  8.50 10.00  6.50
2    0.10  4.00 10.00  6.50
3    0.15  8.50 10.00  6.50
4    0.20  3.50 10.00  6.50
5    0.25  3.50 10.00  6.50
6    0.30  3.00 10.00  6.50
7    0.35  3.00 10.00  6.50
8    0.40  3.00 10.00  6.50
9    0.45  3.50 10.00  0.00
10   0.50  4.50 10.50  0.00
11   0.55  4.50 10.50  0.00
12   0.60  4.50 10.50  0.00
13   0.65  4.50 10.50  0.00
14   0.70  4.50 10.50  0.00
15   0.75  4.50 10.50  0.00
16   0.80  4.50 10.50  0.00
17   0.85  4.50 10.50  0.00
18   0.90  5.00 10.50  0.00
19   0.95  5.50 10.50  0.00
20   1.00  5.50 10.00  0.00
21   1.05  5.50 10.50  0.00
22   1.10  5.50 10.00  0.00
23   1.15  5.50 10.00  0.00
24   1.20  5.50 10.00  0.00
25   1.25  5.50 10.00  0.00
26   1.30  6.00 10.00  0.00
27   1.35  6.00 10.00  0.00
28   1.40  6.50 10.00  0.00
29   1.45  6.50  9.50  0.00
30   1.50  6.50  9.50  0.00
31   1.55  6.00  9.50  0.00
32   1.60  6.00  9.50  0.00
33   1.65  7.00  9.00  0.00
34   1.70  7.00  9.00  0.00
35   1.75  7.00  9.00  0.00
36   1.80  7.00  9.00  0.00
37   1.85  7.00  9.00  0.00
38   1.90  7.50  9.00  0.00
39   1.95  8.00  9.00  0.00
40   2.00  8.50  9.00  0.00
41   2.05  8.50  9.00  0.00
42   2.10  8.50  9.00  0.00
43   2.15  8.50  9.00  0.00
44   2.20  8.50  9.00  0.00
45   2.25  8.50  9.00  0.00
46   2.30  9.00  9.00  0.00
47   2.35  9.50  9.00  0.00
48   2.40  9.50  8.50  0.00
49   2.45  9.50  8.00  0.00
50   2.50 10.00  7.50  0.00
51   2.55 10.00  8.00  0.00
52   2.60 10.50  8.00  0.00
53   2.65 10.50  8.00  0.00
54   2.70 11.00  8.00  0.00
55   2.75 11.50  8.00  0.00
56   2.80 11.50  8.00  0.00
57   2.85 12.00  8.00  0.00
58   2.90 12.00  8.00  0.00
59   2.95 12.00  8.00  0.00
60   3.00 12.00  8.00  0.00
61   3.05 12.00  8.00  0.00
62   3.10 12.00  8.00  0.00
63   3.15 12.00  8.00  0.00
64   3.20 12.00  8.00  0.00
65   3.25 12.00  8.00  0.00
66   3.30 12.00  8.00  0.00
67   3.35 12.00  8.00  0.00
68   3.40 12.00  8.00  0.00
69   3.45 12.00  8.00  0.00
70   3.50 12.00  8.00  0.00
71   3.55 12.00  8.00  0.00
72   3.60 12.00  8.00  0.00
73   3.65 12.00  8.00  0.00
74   3.70 12.00  8.00  0.00
75   3.75 12.00  8.50  0.00
76   3.80 12.00  8.50  0.00
77   3.85 12.00  8.50  0.00
78   3.90 12.00  9.00  0.00
79   3.95 12.00  9.00  0.00
80   4.00 12.00  9.00  0.00

Edit: For the binning per zone part, I want to find the number in which the line segment between start and end points is bisected by the region boundary

Edit (added filter parameters):

Ambulatory Trigger: The minimum number of ambulatory counts since the last stereotypic episode necessary for the current ambulatory movement to be considered an Ambulatory Episode.

Box Size: A user-defined square whose width is measured in beams. The animal must move outside of this box before a movement can be considered ambulatory. Each beam is a 1 inch distance apart. So for maximum distance the "box size" would be a 1 inch square.

Resting Delay: If the subject stays inside the box (see Box Size above) for longer than the resting delay, the movement becomes stereotypic.

jc2525
  • 141
  • 10
  • Binning part - do you want to recreate the whole table with a different duration (say 1.5 seconds between two consecutive rows ) or just want to find distances covered every 1.5 seconds? – R.S. Oct 10 '20 at 18:11
  • I would like it to be variable binning. This sample only shows up to 4 seconds but the whole file could show up to an hour or more. I would like for the user (whoever runs the script) to be able to change the bin type (like, 1 min, 5min, or 10min bins for example) – jc2525 Oct 10 '20 at 23:06

2 Answers2

0

Here is my horrendously long answer. I have never worked with millisecond data and it poses quite a few issues. Plotting a timeseries with millisecond scale axis proved messy too.

Load Data:

    require(raster) #for pointDistance
    require(xts) #For time series based aggregates              
    
    df<-  
      structure(list(indx = 1:80, time = 
                       c(0.05, 0.1, 0.15, 0.2, 0.25, 
                         0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 
                         0.9, 0.95, 1, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3, 1.35, 1.4, 1.45, 
                         1.5, 1.55, 1.6, 1.65, 1.7, 1.75, 1.8, 1.85, 1.9, 1.95, 2, 2.05, 
                         2.1, 2.15, 2.2, 2.25, 2.3, 2.35, 2.4, 2.45, 2.5, 2.55, 2.6, 2.65, 
                         2.7, 2.75, 2.8, 2.85, 2.9, 2.95, 3, 3.05, 3.1, 3.15, 3.2, 3.25, 
                         3.3, 3.35, 3.4, 3.45, 3.5, 3.55, 3.6, 3.65, 3.7, 3.75, 3.8, 3.85, 
                         3.9, 3.95, 4), 
                     
                     x = c(8.5, 4, 8.5, 3.5, 3.5, 3, 3, 3, 3.5, 4.5, 
                           4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 5, 5.5, 5.5, 5.5, 5.5, 5.5, 
                           5.5, 5.5, 6, 6, 6.5, 6.5, 6.5, 6, 6, 7, 7, 7, 7, 7, 7.5, 8, 8.5, 
                           8.5, 8.5, 8.5, 8.5, 8.5, 9, 9.5, 9.5, 9.5, 10, 10, 10.5, 10.5, 
                           11, 11.5, 11.5, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
                           12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12), 
                     y = c(10, 10,10, 10, 10, 10, 10, 10, 10, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5, 
                           10.5, 10.5, 10.5, 10.5, 10, 10.5, 10, 10, 10, 10, 10, 10, 10, 
                           9.5, 9.5, 9.5, 9.5, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 
                           9, 8.5, 8, 7.5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
                           8, 8, 8, 8, 8, 8, 8, 8, 8, 8.5, 8.5, 8.5, 9, 9, 9), 
                     z = c(6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
                     )), 
                class = "data.frame", row.names = c(NA, -80L)) 
    
    
    #remove z
    df$z<- NULL
    head(df)

Find Distance covered per timestep:

    #copy previous location to new column
    df$x_lag<- c(NA,head(df$x,-1))
    df$y_lag<- c(NA,head(df$y,-1))
    
    # Find distance covered in each timestep 
    # lonlat=FALSE for flat surface 
    df$dist<- apply(df, 1 , function(r){pointDistance(c(r[5],r[6] ), c(r[3],r[4] ) , lonlat=FALSE)  })
    head(df)
    sum(df$dist, na.rm = TRUE) #29.44

Find bounding rect for central region:

    ## find bounds for central region 
    box_origin=c(0,0)
    box_size= 17
    
    #center box of size 3
    cent_size= 3 
    cent_start<- ( box_origin[1] + box_size/2 - cent_size/2)
    cent_end<- (box_origin[2] + box_size/2 + cent_size/2)        
    toString(c(cent_start, cent_end))
    
    
    # In the mask below I am totally ignoring any row where the object crosses the boundary
    # It can probably be done much better with a package like sp 
    
    # A mask column to show if the point is in central box
    df$in_center<-(df$x>= cent_start) & (df$x<= cent_end) & (df$y>= cent_start) & (df$y<= cent_end)         
    head(df)        
    table(df$in_center)

Distance Covered :

    #distance covered in center box
    dist_in_center<- sum(df$dist[df$in_center] , na.rm = T)
    
    #add a column with this info
    df$dist_in_center<- 0
    df$dist_in_center[df$in_center]    <- df$dist[df$in_center]
    head(df)
    
    #distance covered in other box
    dist_outof_center<- sum(df$dist[!df$in_center], na.rm = T)
    df$dist_ex_center<- 0
    df$dist_ex_center[!df$in_center]    <- df$dist[!df$in_center]
    
    head(df)
    c(dist_in_center, dist_outof_center)

Plotting : Basically because I don't have unit test

    ######################
    # Plotting 
    ######################
    #plot with polypath
    plot( df[,c(3,4)], type='p',  xlim=c(0, box_size) , ylim= c(0, box_size) , main="Locations", 
          col = ifelse(df$in_center,'red','green') )        
    polypath(df$x, df$y,  rule = "winding" ,border="slateblue4")  #path      
    rect(  cent_start, cent_start, cent_end, cent_end , border = "red" )
    text ( 3,3, toString(floor(c(dist_in_center, dist_outof_center))), col="slateblue4")
    
    #Plot by seconds
    plot( df[,c(3,4)], type='p',  xlim=c(0, box_size) , ylim= c(0, box_size) , main="Locations", 
          col = floor(df$time) )
    rect(  cent_start, cent_start, cent_end, cent_end , border = "red" )
    text ( 3,3, toString(floor(c(dist_in_center, dist_outof_center))), col="slateblue4")

Locations and Bounding Box

Location by Seconds

Time Series:

    ##########################
    #Convert to Time Series and aggregate by time periods 
    ##########################
    
    
    require(xts)
    
    #make R print milliseconds 
    #https://stackoverflow.com/questions/9778125/r-xts-millisecond-index
    opts <- options(digits.secs = 3)
    getOption("digits.secs")
    # This does not work for timestamps ending in :00 
    #as.POSIXct(as.character(df$time), format = "%S.%OS")  
    #as.POSIXct(as.character(df$time))
    
    # This works for all 
    as.POSIXct(format(df$time, nsmall = 3), format = "%S.%OS") 
    
    colnames(df)
    
    dist_xts<- xts(df[-1,c("dist", "dist_in_center","dist_ex_center")], order.by = as.POSIXct(format(df$time[-1], nsmall = 3), format = "%S.%OS") )
    
    
    # Sum per minute 
    xts::period.sum(dist_xts$dist, xts::endpoints(dist_xts$dist, on="minutes"))
    
    # Sum per second 
    xts::period.sum(dist_xts$dist, xts::endpoints(dist_xts$dist, on="seconds"))
    xts::period.sum(dist_xts$dist_in_center, xts::endpoints(dist_xts$dist_in_center, on="seconds"))
    xts::period.sum(dist_xts$dist_ex_center, xts::endpoints(dist_xts$dist_in_center, on="seconds"))
    
    tail(df)

Without a Time Series :

    # ######################
    # # Aggregating distance over 3 timesteps , without time series 
    # ######################
    # 
    # 
    # sum3_T<- c(rep(NA,2),  rollapply(df$dist*as.numeric(df$in_center), 3, sum))
    # sum3_F<- c(rep(NA,2),  rollapply(df$dist*as.numeric(!df$in_center), 3, sum))
    # sum3_T
    # sum3_F
    # 
    # # distance covered in last 3 timesteps
    # # df$sum3_overall<- sum3_T+sum3_F
    # df$sum3_overall<- c(rep(NA,2),  rollapply(df$dist, 3, sum))
    # # distance covered in last 3 timesteps in current zone only
    # df$sum3_zoned<- ifelse(df$in_center, sum3_T , sum3_F)
    # df$sum3_outofzone<- df$sum3_overall - df$sum3_zoned
    # 
    # head(df)
R.S.
  • 2,093
  • 14
  • 29
  • Thank you! That was a lot of help. It's gonna take me some time to read through it. Why did you ignore rows that cross the boundary? – jc2525 Oct 11 '20 at 16:44
  • I ignored crossovers because then we will have to make some decisions , which only you can make. For example- 1. Do we assing the distance covered to starting region or the ending region? (solution is simple, just use half the masking logic used above ) OR 2. Do we split it by half and assign to both regions , OR 3. Do we actually find the ratio in which the line segment between start and end points is bisected by the region boundary. (there must be a package for that) – R.S. Oct 11 '20 at 20:00
  • I understand. I am looking more for scenario 3 unfortunately. – jc2525 Oct 11 '20 at 20:23
  • it sounds interesting , I'll look into that. Maybe you should add this requirement to the question. It might attract the attention of someone who knows the solution already – R.S. Oct 11 '20 at 20:40
  • I added 3 parameters to my question to help. The company software that supplies the output would not give me the algoritm but did tell me it was based off of the parameters. Here is the manual...https://www.med-associates.com/wp-content/uploads/2017/01/DOC-038-R3.2-SOF-811-ACTIVITY-MONITOR.pdf – jc2525 Oct 16 '20 at 19:43
  • @jc2525 I think I did solve the ratio problem that day but then got busy. Wasn't very elegant. Was looking for either a straightforward spatial library that could handle everything, or a hands-on linear solver. Will post what I have got, later in the day. – R.S. Oct 18 '20 at 02:09
  • Thank you for keeping up with this! I really appreciate it :) – jc2525 Oct 19 '20 at 23:53
  • Posting a new Answer. I haven't looked into the conditional ambulation , but I think it will be relatively easy for thereon. – R.S. Oct 20 '20 at 08:32
0

Adding a new answer with a function to find exact ratios in which a path is segmented by the polygon

Summary:

  • Store lines and rectangle as SpatialLines SpatialPolygons
  • Create a function to find ratios in which a polygon intersects a line
  • lapply the function to all consecutive pairs of points
  • reduce the list into a dataframe
  • convert to xts and aggregate by desired timesteps

Load Libraries

require(xts) #For time series based aggregates              
require(sp) #for spatial data structures
require(raster)#for making computations on spatial data

Load Data

df<-  
  structure(list(indx = 1:80, time = 
                   c(0.05, 0.1, 0.15, 0.2, 0.25, 
                     0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 
                     0.9, 0.95, 1, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3, 1.35, 1.4, 1.45, 
                     1.5, 1.55, 1.6, 1.65, 1.7, 1.75, 1.8, 1.85, 1.9, 1.95, 2, 2.05, 
                     2.1, 2.15, 2.2, 2.25, 2.3, 2.35, 2.4, 2.45, 2.5, 2.55, 2.6, 2.65, 
                     2.7, 2.75, 2.8, 2.85, 2.9, 2.95, 3, 3.05, 3.1, 3.15, 3.2, 3.25, 
                     3.3, 3.35, 3.4, 3.45, 3.5, 3.55, 3.6, 3.65, 3.7, 3.75, 3.8, 3.85, 
                     3.9, 3.95, 4), 
                 
                 x = c(8.5, 4, 8.5, 3.5, 3.5, 3, 3, 3, 3.5, 4.5, 
                       4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 5, 5.5, 5.5, 5.5, 5.5, 5.5, 
                       5.5, 5.5, 6, 6, 6.5, 6.5, 6.5, 6, 6, 7, 7, 7, 7, 7, 7.5, 8, 8.5, 
                       8.5, 8.5, 8.5, 8.5, 8.5, 9, 9.5, 9.5, 9.5, 10, 10, 10.5, 10.5, 
                       11, 11.5, 11.5, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
                       12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12), 
                 y = c(10, 10,10, 10, 10, 10, 10, 10, 10, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5, 
                       10.5, 10.5, 10.5, 10.5, 10, 10.5, 10, 10, 10, 10, 10, 10, 10, 
                       9.5, 9.5, 9.5, 9.5, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 
                       9, 8.5, 8, 7.5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
                       8, 8, 8, 8, 8, 8, 8, 8, 8, 8.5, 8.5, 8.5, 9, 9, 9), 
                 z = c(6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
                 )), 
            class = "data.frame", row.names = c(NA, -80L)) 
#remove z
df$z<- NULL
head(df)
tail(df)
dim(df)

Path as Spatial Data: (Not needed, but good to plot)

line.path<- sp::Line( df[,c("x","y")] ) 
#convert to SpatialLines which give more options
spatlines.path= sp::SpatialLines( list(Lines(list(line.path), ID="path")) )
summary(spatlines.path)
plot(spatlines.path, main="path", xlim=c(1,13), ylim=c(5,12))

Create A Polygon that acts as our central region of interest: I am using an arbitrary rectangle, You can compute its exact coordinates,just like in my previous answer. We will use this polygon throughout this code.

This should also work for other shapes

######### POLYGON ########################
#let these be the coordinates of our critical/central  region :
#( compute bounds  as per your requirement)

rect.df<- data.frame(x=c(6,10,10,6), y=c( 7,7,11,11)) 
poly.rect<- sp::Polygon(coords = rect.df)
# you can find which points lie on this polygon . Not using this 
# out=0, in=1, onboundary=2
sp::point.in.polygon( point.x = df$x,point.y = df$y, 
                      pol.x = poly.rect@coords[,1], pol.y = poly.rect@coords[,2] )

#Convert polygon to  spatial polygons for more options. 
#?`SpatialPolygons-class`
spat.poly<- sp::SpatialPolygons( list(Polygons(list(poly.rect), ID="rectpoly")) )
#to add to previous plot
lines(poly.rect,col="red")

A Function for Line Polygon Intersection

This gives us the exact ratios/segments in which boundaries of a polygon divides a line.

#@param pts= representation of two points that constitute the line
#@param poly= Spatial Polygon
#@returns c( total_length, inside_len, outside_len )
line.polygen.intersects<- function ( pts = df[3:4,c("x","y")] , poly ){
  
  # print(pts)
  # if both points are same
  if (all(( pts[1,]-pts[2,])==0)){ return(c(0,0,0))}
  
  total.length<- LineLength(as.matrix(pts))
  
  line1<- sp::Line(coords = pts)
  spat.line1<- sp::SpatialLines( list(Lines(list(line1), ID="lineseg")) )
  
  #crop line as per polygon boundaries
  crop.line.poly<- raster::crop(spat.line1,  poly)
  #totally outside , crop is null
  if(is.null(crop.line.poly)){return(c(total.length,0,total.length ))}
  #if crop is just one point, coords are vector
  if(!is.matrix( coordinates(crop.line.poly)[[1]][[1]] ) ){
                       return(c(total.length,0,total.length ))}
  
  #print(coordinates(crop.line.poly))
  #print(crop.line.poly)
  in.length<- LineLength( coordinates(crop.line.poly)[[1]][[1]],longlat = FALSE)
  
  plot(NULL, main="plot", xlim=c(0,15), ylim=c(0,15))
  points(poly@polygons[[1]]@Polygons[[1]]@coords, type = "l", col="red")
  points(pts, type="l")
  
  return( c( total.length, in.length, (total.length-in.length) ))
}

#spat.poly
line.polygen.intersects(poly=spat.poly)
line.polygen.intersects(pts = matrix(c(0,8,0,8),2,2,byrow =F),poly= spat.poly)
#0 length
line.polygen.intersects(pts = matrix(c(0,0, 8,8),2,2,byrow = F),poly= spat.poly)
#line totally outside polygon .. this is problematic, else are ok
line.polygen.intersects(pts = matrix(c(0,0,0,8),2,2, byrow = F),poly= spat.poly)
#line totally inside polygon
line.polygen.intersects(pts = matrix(c(7,9,8,10),2,2, byrow=F),poly= spat.poly)
#line intersecting at two points 
line.polygen.intersects(pts = matrix(c(0,10,0,12),2,2,byrow=F),poly= spat.poly)
#line overlapping with one side 
line.polygen.intersects(pts = matrix(c(0,12,7,7),2,2,byrow=F),poly= spat.poly)
# when both points are bang on a side of poly
line.polygen.intersects(pts = matrix(c(6,10,10,8),2,2,byrow=F),poly=spat.poly)
# when exactly 1 point on edge and rest is outside 
line.polygen.intersects(pts = df[25:26,c("x","y")], poly= spat.poly)

Compute Intersections for all data points

seglengths<- lapply( 2:80,  function(i){ print(i); 
               line.polygen.intersects(df[(i-1):i, c("x","y")], poly = spat.poly)})
#convert to dataframe
df_distances<- data.frame(Reduce( rbind , seglengths) )
colnames(df_distances)<- c("len.total","len.in","len.out")
df_distances<- rbind( data.frame(len.total=0,len.in=0,len.out=0), df_distances)
rownames(df_distances)<- row.names(df) 
View(df_distances)

apply(df_distances,2,sum)
#merge with original data
df_distances<- cbind(df[,c("time","x","y")], df_distances)
df_distances

Convert to TimeSeries

###########################################
# convert to Time series and aggregate timewise
###########################################
require(xts)

opts <- options(digits.secs = 6)
getOption("digits.secs")
colnames(df_distances)

dist_xts<-  xts(df_distances, as.POSIXct(df$time, origin="2000-01-01", tz="UTC") )

View(dist_xts)

Aggregate by Time

###########################################
# Aggregate by Time Periods
###########################################


# Sum byperiod
xts::period.sum(dist_xts$len.total, xts::endpoints(dist_xts$len.total, on="minutes"))
xts::period.sum(dist_xts$len.total, xts::endpoints(dist_xts$len.total, on="seconds"))

# aggregate  with truncation handles the index better :  
aggregate(dist_xts$len.total, as.POSIXct(trunc(time(dist_xts), "sec")), sum)
aggregate(dist_xts$len.total, as.POSIXct(trunc(time(dist_xts), "sec")), mean)
?aggregate
?trunc.POSIXt
aggregate(dist_xts$len.total,by = as.POSIXct(trunc(time(dist_xts), units="secs")), sum)
aggregate(dist_xts$len.total,by = as.POSIXct(trunc(time(dist_xts), "min")), sum)
aggregate(dist_xts[,-1], as.POSIXct(trunc(time(dist_xts), "sec")), sum)

Output:

                        x     y len.total   len.in  len.out
2000-01-01 00:00:00  87.0 195.0 17.118034 7.500000 9.618034
2000-01-01 00:00:01 127.0 191.5  5.618034 3.618034 2.000000
2000-01-01 00:00:02 199.5 168.0  5.707107 3.707107 2.000000
2000-01-01 00:00:03 240.0 163.5  1.000000 0.000000 1.000000
2000-01-01 00:00:04  12.0   9.0  0.000000 0.000000 0.000000

# Some fractional seconds appear rounded off oddly. but that doesn't affect computation
# This discussion might be useful
# https://stackoverflow.com/questions/7726034/how-r-formats-posixct-with-fractional-seconds
R.S.
  • 2,093
  • 14
  • 29