Counting How Many Times a Point Passes Through a Grid of Coordinates

Hi Everyone,
I am currently trying to figure out how to count the number of times an object passes through each grid point defined by a set of coordinates. My grid is defined by x-coordinates and y-coordinates as:

longs = long[1]:Δ:long[2]
lats = lat[2]:-Δ:lat[1]
xx = longs' .* ones(length(lats))
yy = ones(length(longs))' .* (lats)


long = [-95.0, -75.0]
lat = [50.0, 62.5]
Δ = 0.005

What want to do is create a grid of x and y coordinates along those latitude and longitudes with a spacing of 0.005deg (which I believe I have successfully created). I now have traffic that transits along a path with a certain lat/long and I want to count the number of times each of the traffics positions fall within each grid cell but I am not sure how to do this. Would anyone be able to provide some guidance on how I can go about doing this or if there is a package that can do this? I will be counting the number of times traffic passes through each cell for ~500k different transits.

An example of a small transit path would be:

transit = Dict("longs"=>[-87.84, -85.03, -82.69, -81.54, -81.35], "lats"=>[61.57, 59.83, 56.51, 54.90, 54.29])

The following example was produced for fun only and it is certainly not the most efficient but I believe it works within a given error tolerance (controlled by the number of linearly interpolated parametric spline points N4). It assumes that each trajectory does not do repeated passes over the same grid cells.

Then one just needs to set a counter with zeros over the lat-long grid and update it with the unique grid cell intersections found for each trajectory. In the code below this would mean increasing the counter by one at each grid entry given by the integer indices: xyi[1,i], xyi[2,i] .

PS: updated code with a more logical definition of parameter N4

using Dierckx, Plots; gr(legend=false)

long = [-90.0, -80.0]
lat = [53.0, 63.0]
Δ = 0.5

longs = long[1]:Δ:long[2]
lats = lat[1]:Δ:lat[2]

transit = Dict("longs"=>[-87.84, -85.03, -82.69, -81.54, -81.35], "lats"=>[61.57, 59.83, 56.51, 54.90, 54.29])
tlon = transit["longs"]
tlat = transit["lats"]
scatter(tlon, tlat, mc=:blue, grid=false, ms=3)
vline!(longs, lc=:lightgrey,lw=0.3); hline!(lats, lc=:lightgrey,lw=0.3)

N = length(tlon);
# Increase N4 for better accuracy if needed
N4 = 50*Int( (diff([extrema(tlon)...])[1]÷Δ + 1)*(diff([extrema(tlat)...])[1]÷Δ + 1) )
t = LinRange(0, 1, N)
spl = ParametricSpline(t, [tlon tlat]', k=1)  # linear splines
tfine = LinRange(0, 1, N4)
xys = evaluate(spl,tfine);
plot!(xys[1,:], xys[2,:], lc=:red, lw=1)

mid_xy = fill(Int(0), 2, N4)
@. mid_xy[1,:] = (xys[1,:] - long[1])÷Δ + 1
@. mid_xy[2,:] = (xys[2,:] -  lat[1])÷Δ + 1

xyi = unique(mid_xy,dims=2)
Nint = size(xyi,2)

# Display intersections found
println("Sampling trajectory at $N4 points found $Nint grid cell intersections")
rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
for i in 1:Nint
    display(plot!(rectangle(Δ, Δ, longs[xyi[1,i]], lats[xyi[2,i]]), opacity=.1, color=:blue))
plot!(title = "Number of grid cell intersections found ~ $Nint")



Here is another implementation (that I made for a similar problem a long time ago). Same caveats apply though: could be more efficient, could be neater! It doesn’t use any packages though, which may be desirable.

function segment_line_2d(p0, p1, d;
                         include_startpoint = false,
                         include_endpoint   = false)

    x0, y0 = p0
    x1, y1 = p1

    if d isa Number
        d = (d, d)

    # Coordinates at which a grid border is crossed:
    x_crossings = border_crossings(x0, x1, d[1])
    y_crossings = border_crossings(y0, y1, d[2])

    # Complementary coordinates for each of the crossings:
    x_complements_to_y = complements(y0, y1, x0, x1, y_crossings)
    y_complements_to_x = complements(x0, x1, y0, y1, x_crossings)

    all_xs = Iterators.flatten((x_crossings, x_complements_to_y))
    all_ys = Iterators.flatten((y_complements_to_x, y_crossings))

    # Sort the points in order from (x0,y0) to (x1,y1)
    points = sort!(collect(zip(all_xs, all_ys)), rev = x0 > x1)

    include_startpoint && pushfirst!(points, Tuple(p0))
    include_endpoint && push!(points, Tuple(p1))

    return unique!(points)

# return a range including all of the intersection values in 
# the exclusive range x0:x_spacing:x1 (left-to-right vs 
# right-to-left sensitive)
function border_crossings(x0, x1, x_spacing)
    x0_pixel = floor(Int64, x0/x_spacing)
    x1_pixel = floor(Int64, x1/x_spacing)

    n_xs = abs(x1_pixel - x0_pixel)

    # left-to-right vs right-to-left
    if x1_pixel > x0_pixel
        return (x0_pixel .+ (1:n_xs)) .* x_spacing
        return (x0_pixel .+ (1 .- (1:n_xs))) .* x_spacing

# use 2-point line formula to return corresponding 
# y values for known x-grid-crossings
function complements(x0, x1, y0, y1, x_crossings)
    m = (y1 - y0) / (x1 - x0)
    return m .* (x_crossings .- x0) .+ y0

transit = Dict("longs"=>[-87.84, -85.03, -82.69, -81.54, -81.35],
               "lats"=>[61.57, 59.83, 56.51, 54.90, 54.29])

pts = collect(zip(transit["longs"], transit["lats"]))

inters = NTuple{2, Float64}[]
for i in 1:length(pts)-1
    append!(inters, segment_line_2d(pts[i], pts[i+1], Δ))

using Plots
plot(pts, label = "waypoints")
vline!(-88:0.5:-81, lc=:lightgrey,lw=0.3, label = nothing)
hline!(54:0.5:62, lc=:lightgrey,lw=0.3, label = nothing)
scatter!(inters, label = "intersections")


That said, I find the wording of your question a bit confusing, so I’m wondering if this or the previous reply answer it. Particularly:

What does the section I’ve bolded mean exactly?

1 Like

I think both of these answers should provide a start at least for what I want to do. Essentially what I am doing is looping through all transits and counting the number of transits that pass through each grid cell, that is what I tried to convey in the bolded statement.

I ended up converting the following code from MATLAB to Julia and it works great:

I was having some problems with interpolation as the points are necessarily increasing in x dimension but this solution does the job and was quite easy to convert to Julia.