Check if a point is near the country boundary

General talk
1

Hi all and happy new year!
We hope 2021 will be better than 2020, or at least not badder… :smiley:

My problem: I have a list of points (IFR reporting point, from Eurocontrol) in Europe and other countries, but I’m just interesting on the points in Europe.
I’d like to filter them and just get the ones near (let’s say 20 Km) the country boundaries, since these are the point that VFR-pilots uses to report entering a new country.
I imported the OSM data for Europe in my PostgreSQL database and I now I have the data to display the country boundaries.
But I really don’t have an idea how to build an SQL query to check if a give point (latitude + longitude) is 20 Km from the line of the country boundary and return TRUE/FALSE and the affected countries.

Can someone help me?

Thanks a lot
Luca

2

Hi again,

maybe I found something, but it does not give right data…

I see the function ST_Distance. This should do what I need, but if I try to calculate the distance from a point near the boundary between Germany and Czech republic, I get wrong values:

SELECT
  name,
  admin_level,
  ST_Distance(ST_Transform('SRID=4326;POINT(50.8379138888889 14.0007583333333)'::geometry, 3857), way) * cosd(42.3521) as dist
FROM planet_osm_roads
WHERE boundary = 'administrative'

  AND admin_level = '2'
  AND osm_id < 0
ORDER BY admin_level DESC;


name     | admin_level |       dist

-------------±------------±-----------------
Česko | 2 | 4812054.73000158
Deutschland | 2 | 4812054.73000158
Česko | 2 | 4812014.2751998
Deutschland | 2 | 4812014.2751998
Česko | 2 | 4801322.76686319
Deutschland | 2 | 4801322.76686319
Česko | 2 | 4798301.61335419
Deutschland | 2 | 4793808.13783476
Česko | 2 | 4786576.73402148
Deutschland | 2 | 4785266.87451821
Česko | 2 | 4775493.77594454
Deutschland | 2 | 4775493.77594454
Česko | 2 | 4746495.15621187
Deutschland | 2 | 4746495.15621187
Česko | 2 | 4775724.71111992
Deutschland | 2 | 4767387.92338493
Polska | 2 | 4748955.14857684
Polska | 2 | 4759696.4876026
Deutschland | 2 | 4759696.4876026
Deutschland | 2 | 4849836.9269982
Polska | 2 | 4849836.9269982
Polska | 2 | 4797231.65721602
Deutschland | 2 | 4797231.65721602
Polska | 2 | 4890476.28917314
Deutschland | 2 | 4890476.28917314
Polska | 2 | 4938293.09438962
Deutschland | 2 | 4938293.09438962
Deutschland | 2 | 5002527.32842157
Polska | 2 | 5002527.32842157
(29 Zeilen)

can someone explain me what I’m doing wrong?

Thanks
Luca

3

What you want to do is a) buffer the boundary by a given amount (say 100 m) using ST_BUFFER & b) then do an ST_INTERSECTS between your roads & the boundary & c) only do distances on the roads which intersect. You may find it more efficient to save the buffered data as a working table which can be indexed for efficiency, particularly as you dont have things in the projection you are using.

Note this approach works symmetrically both sides of a boundary, it may be more complex if you want things only on one side of the boundary.

The basic idea is that you want to reduce total query data fairly quickly & enable use of indexes. Personally I’d use ETRS-89 as the projection as it is better than Popular Mercator for Europe.

4

Hi and thank you for your answer!
I’m feeling really dumb, but I can’t understand what you mean…
First, I don’t have “roads”, but “points”. Just point, identified by latitude and longitude.
Then I can’t understand what you mean with the intersection…

Maybe could you give me an example as SQL to explain what you mean?

Thank you very much!
Luca

5

Hi again,

maybe I undestand what you mean… something like that?

SELECT name, ST_Intersects(ST_Buffer(ST_GeomFromText('POINT(50.8379138888889 14.0007583333333)'), 20000, 'quad_segs=8'), way)
FROM planet_osm_roads
WHERE boundary = 'administrative'
  AND admin_level = '2'
  AND osm_id < 0;

Unfortunately it does not work as expected, since all data return false:

name     | st_intersects

-------------±--------------
Deutschland | f
Česko | f
Česko | f
Deutschland | f
Česko | f
Deutschland | f
Česko | f
Deutschland | f
Česko | f
Deutschland | f
Česko | f
Deutschland | f
Česko | f
Deutschland | f
Česko | f
Deutschland | f
Polska | f
Polska | f
Deutschland | f
Polska | f
Deutschland | f
Polska | f
Deutschland | f
Deutschland | f
Polska | f
Polska | f
Deutschland | f
Deutschland | f
Polska | f
(29 Zeilen)

What do I do wrong?

Thanks
Luca

6

No something like this (not checked), assuming planet_osm_roads is in popular Mercator & you have points in a table in, say, WGS84

SELECT ifr_points.*, planet_osm_roads.“name”, st_distance(st_transform(geom,3857), way) distance
FROM planet_osm_roads, ifr_points
WHERE st_intersects(st_transform(ifr_points.geom,3857), st_buffer(way,20000))
and boundary = ‘administrative’
AND admin_level = ‘2’
AND osm_id < 0;

7

I got it with this query:

SELECT name, MIN(dist) FROM
(
  SELECT
    name,
    ST_Distance('SRID=4326;POINT(14.0007583333333 50.8379138888889)'::geography, st_transform(way,4326)::geography) as dist
  FROM planet_osm_roads
  WHERE boundary = 'administrative'
    AND admin_level = '2'
    AND osm_id < 0
  ORDER BY dist ASC
) a
WHERE dist < 20000
GROUP BY name;

My problem was that “POINT” need first the longitude and then the latitude… :stuck_out_tongue:

Bye
Luca

8

Hi Luca,

For a single point that’s fine. I thought you had lots!

Jerry

9

Hi Jerry,

correct! I have a list with many hundert of points. I’ll insert them in a table and JOIN it with planet_osm_roads, of course.
What I needed was just to check if what I want is possible and how to do that…

Now the server import the last data from OSM. As far it is done, I’ll insert the data in a table and get all boundary’s reporting points.

10

Hi Luca,

You may find performance different with lots of points, but YMMV.

Jerry