Stream: helpdesk (published)

Topic: ✔ How to check if position `(lon, lat)` is within an area?

Sascha Mann (Mar 28 2025 at 18:39):

Hey, this is probably a stupid question but I'm completely new to working with map data and all the various geo packages are a bit confusing.

I have a shapefile that contains several districts and want to check which district a point is in given by latitude and longitude of the point.

I found a thread on Discourse that has several options but I'm at a loss which one to use. I've tried Júlio Hoffimann's method but it didn't match any points (which is probably me making a mistake somewhere):

geotable = GeoIO.load("Kommunalwahlbezirk_3/Kommunalwahlbezirk.shp")
ps = Point.(points)
gs = domain(geotable)
[p ∈ g for p in ps, g in gs]

with

julia> ps = Point.([(lon, lat)])
1-element Vector{Point{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}:
 Point(x: 6.958137997831819 m, y: 50.941303500000004 m)
 ...

and

julia> gs
45 GeometrySet
├─ Multi(2×PolyArea)
├─ Multi(1×PolyArea)
├─ Multi(1×PolyArea)
├─ Multi(1×PolyArea)
├─ Multi(1×PolyArea)
⋮
├─ Multi(1×PolyArea)
├─ Multi(1×PolyArea)
├─ Multi(1×PolyArea)
├─ Multi(1×PolyArea)
└─ Multi(1×PolyArea)

However the resulting matrix has no matches, even though that point is definitely within the area of one of the districts.

Any ideas on what I'm doing wrong with that specific method or generally what the best idea to approach this problem is? It doesn't have to be efficient, I only have to do this once for ~2.3k points.

Daniel González (Mar 28 2025 at 19:20):

You probably need to define the Point in the same projection, I'll try to see if I can find out how to construct it

Sascha Mann (Mar 28 2025 at 19:22):

Daniel González said:

You probably need to define the Point in the same projection, I'll try to see if I can find out how to construct it

I literally just arrived at that thought when you wrote it :sweat_smile: the shapefile coordinates are significantly off from the coordinates I have. The coordinate system for the shapefile is ETRS 1989 UTM Zone 32N

Sascha Mann (Mar 28 2025 at 19:23):

The coordinates of my points are what I receive from the OpenStreetMaps API (which appears to be WGS 1984)

Daniel González (Mar 28 2025 at 19:31):

I think this might be what you need

using CoordRefSystems, Unitful

m = Unitful.m
° = Unitful.°

k₀ = 0.9996
latₒ = 0.0°
lonₒ = 9.0°
xₒ = 500000.0m
yₒ = 0.0m
S = CoordRefSystems.Shift(; lonₒ, xₒ, yₒ)
proj = TransverseMercator{k₀,latₒ,ETRFLatest,S}

point_projected = convert(proj, LatLon(lat, lon)) # use your trial values here

Daniel González (Mar 28 2025 at 19:35):

Alternatively

proj = utm(:north, 32; datum=ETRFLatest)

Sascha Mann (Mar 28 2025 at 19:49):

That appears to work, thank you!

Notification Bot (Mar 28 2025 at 20:04):

Sascha Mann has marked this topic as resolved.

Júlio Hoffimann (Mar 31 2025 at 13:22):

@Sascha Mann try to use the Proj transform instead as described in Chapter 6 Map projections: https://juliaearth.github.io/geospatial-data-science-with-julia/06-projections.html

Júlio Hoffimann (Mar 31 2025 at 13:23):

You can send the geotable directly to it as in geotable |> Proj(utm(:north, 32, datum=ETRFLatest))

Júlio Hoffimann (Mar 31 2025 at 13:25):

Thank you @Daniel González for helping with it :100:

Last updated: Apr 25 2025 at 04:43 UTC