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

**URL:** <https://discourse.julialang.org/t/counting-how-many-times-a-point-passes-through-a-grid-of-coordinates/57229>\
**Category:** General Usage\
**Tags:** question\
**Created:** [March 15, 2021, 7:52pm UTC](https://discourse.julialang.org/t/counting-how-many-times-a-point-passes-through-a-grid-of-coordinates/57229 "2021-03-15T19:52:09Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![PJoynt](https://avatars.discourse-cdn.com/v4/letter/p/f4b2a3/32.png) [@PJoynt](https://discourse.julialang.org/u/PJoynt)\
**Post date:** [March 15, 2021, 7:52pm UTC](https://discourse.julialang.org/t/counting-how-many-times-a-point-passes-through-a-grid-of-coordinates/57229/1 "2021-03-15T19:52:09Z")

</div>

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:

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

```

Where:

```julia
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:

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

```

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [March 15, 2021, 10:53pm UTC](https://discourse.julialang.org/t/counting-how-many-times-a-point-passes-through-a-grid-of-coordinates/57229/2 "2021-03-15T22:53:44Z")

</div>

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`_

```julia
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))
end
plot!(title = "Number of grid cell intersections found ~ $Nint")

```

![Trajectory_intersections_on_grid_cells](https://global.discourse-cdn.com/julialang/original/3X/4/9/490e7d4b7ca650ff3f61904292da125bf1251d6f.png)

---

<div class="post-metadata">

**Author:** ![tomerarnon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomerarnon/32/3170_2.png) [@tomerarnon](https://discourse.julialang.org/u/tomerarnon)\
**Post date:** [March 16, 2021, 12:58pm UTC](https://discourse.julialang.org/t/counting-how-many-times-a-point-passes-through-a-grid-of-coordinates/57229/3 "2021-03-16T12:58:27Z")

</div>

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.

```julia

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)
    end

    # 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)
end

# 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
    else
        return (x0_pixel .+ (1 .- (1:n_xs))) .* x_spacing
    end
end

# 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
end

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], Δ))
end

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")

```

![discourse_segment2d](https://global.discourse-cdn.com/julialang/original/3X/a/b/ab397eee6be7d202ce8494fc349f449ea2ff088a.png)

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:

> [@PJoynt](#):
>
> 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.

What does the section I’ve bolded mean exactly?

---

<div class="post-metadata">

**Author:** ![PJoynt](https://avatars.discourse-cdn.com/v4/letter/p/f4b2a3/32.png) [@PJoynt](https://discourse.julialang.org/u/PJoynt)\
**Post date:** [March 16, 2021, 2:54pm UTC](https://discourse.julialang.org/t/counting-how-many-times-a-point-passes-through-a-grid-of-coordinates/57229/4 "2021-03-16T14:54:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![PJoynt](https://avatars.discourse-cdn.com/v4/letter/p/f4b2a3/32.png) [@PJoynt](https://discourse.julialang.org/u/PJoynt)\
**Post date:** [March 18, 2021, 6:49pm UTC](https://discourse.julialang.org/t/counting-how-many-times-a-point-passes-through-a-grid-of-coordinates/57229/5 "2021-03-18T18:49:45Z")

</div>

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

> **[Fast and Robust Curve Intersections](https://www.mathworks.com/matlabcentral/fileexchange/11837-fast-and-robust-curve-intersections)**
>
> Computes intersection points of two curves.

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.