One possible solution is to use the `Aggregate`

transform from GeoStats.jl:

```
Aggregate(domain, var₁ => agg₁, var₂ => agg₂, ..., varₙ => aggₙ)
Aggregate([g₁, g₂, ..., gₙ], var₁ => agg₁, var₂ => agg₂, ..., varₙ => aggₙ)
Aggregate variables var₁, var₂, ..., varₙ over geospatial domain using aggregation functions agg₁, agg₂, ..., aggₙ. Alternatively, aggregate variables over geometries g₁,
g₂, ..., gₙ. Default aggregation function is mean for continuous variables and first otherwise.
Examples
≡≡≡≡≡≡≡≡
Aggregate(domain, 1 => last, 2 => maximum)
Aggregate(domain, :a => first, :b => minimum)
Aggregate(domain, "a" => last, "b" => maximum)
Aggregate(geoms, "a" => last, "b" => maximum)
```

Pass your own aggregation function `aggfun`

to ignore `NaN`

values:

```
img = GeoIO.load("foo.jpg")
grid = CartesianGrid(M, N)
img |> Aggregate(grid, "color" => aggfun)
```

Of course this solution is mostly useful if you are already taking advantage of the full stack we are building for geospatial data science. Others can suggest simpler solutions with smaller packages.

In case you want to learn more: Geospatial Data Science with Julia