Checking that a position is inside a AREA works in EPSG:25382 projection, but not in CRS84

Specific Domains Geo
1

I have first to say, that I’m a newbie when it comes to geo-related stuff, so please excuse me if I ask or formulate something stupid here, but I’ve been biting my teeth out on this problem for hours.

First a little context: I have a table with coordinates in EPSG:3035 and two different shapefiles, the first one is for all of Germany, and if a coordinate is in Berlin, I want the exact district. The shapefiles have different coordinate systems. After some back and forth I was in favor of the first one (germany), but the same methodology produces an error with the berlin shapefile. Here is the MWE with additional comments.

import GeoIO: LatLon, Point, load

# First download the shapefiles used in the MWE and unzip it
# (unzip must be installed on the system)
zipfile = download("https://box.fu-berlin.de/s/wk9ixTJGoHgWSA5/download/shapefiles.zip")
run(`unzip $(zipfile) -d /tmp`)

# the first shapefile uses EPSG:25832 projection and extract the area
# of the first row of the geotable
germany = load("/tmp/kfz250.utm32s/KFZ250.shp")
g_area = germany[1, :geometry]
# we create a point in this projection and can check if this point is
# inside the area
@assert (Point(3.2e5, 5.4e6) in g_area) == false
@assert (Point(3.7302e5, 5.612e6) in g_area) == true

# we try to do the same with another shapefile, which is in CRS84
berlin = load("/tmp/bezirksgrenzen_berlin.shp.zip")
b_area = berlin[1, :geometry]

# when we create a point in CRS84 using LatLon and check if this is in the area
# this will raise this error:
# ERROR: MethodError: no method matching iterate(::Meshes.PolyArea{Meshes.🌐, CoordRefSystems.GeodeticLatLon{…}, Meshes.Ring{…}, Vector{…}})
LatLon(52.5, 13) in b_area

So my question is: Why is the in function call working for the first shapefile but not for the second? And how can I fix this?

Thank you,
Steffen

1 Like
2

Hi @sfuerst , welcome to the forum.

We distinguish between the concept of Point and the concept of coords of the Point in our software stack. Check the outputs in your example:

julia> Point(3.2e5, 5.4e6)
Point with Cartesian{NoDatum} coordinates
├─ x: 320000.0 m
└─ y: 5.4e6 m

julia> LatLon(52.5, 13)
GeodeticLatLon{WGS84Latest} coordinates
├─ lat: 52.5°
└─ lon: 13.0°

The first object is a Point with Cartesian coordinates without a specified datum. The second object are GeodeticLatLon coordinates with WGS84Latest datum.

You can construct Point with coordinates of any type:

julia> Point(LatLon(52.5, 13))
Point with GeodeticLatLon{WGS84Latest} coordinates
├─ lat: 52.5°
└─ lon: 13.0°

The point in geom predicate does not work if you only pass coordinates. It expects a Point.

Also, when you type Point(3.2e5, 5.4e6) in g_area you are assuming that g_area is in Projected or Cartesian coordinates. If your shapefiles come in LatLon coordinates, it would be safer to create a point in the same coordinate reference system. You can extract the crs(g_area) or the crs(germany) in your example to construct coordinates of the same type.

Please let me know if you need further assistance. Happy to help you with these concepts.

A good resource to learn this in more depth is the GDSJL book:

2 Likes