5

I have a point A (EPSG 4326) in Berlin, Germany:

{
   "type": "Point",
   "coordinates": [
     13.394994735717773,
     52.514849775310154
   ]
}

I want to create point B that is 500m to the east of point A and point C that is 500m to the North.

For this, I have to project point A to a projection that preserves distance.

Can you recommend a projection that works for this operation in entire Germany?

I am working with GeoPandas and Shapely for this project.

  • 2
    What are the requierements for the accuracy? For 500m it might be negligible to use the UTM Zone 33 as the deviations around 13.4° are around 25cm per kilometer for this projection. – Matte Apr 23 '17 at 9:04
  • 25cm per kilometer is totally fine. Minimum requirement would roughly be 1m per km. – ustroetz Apr 23 '17 at 9:38
6

You don't need an equidistant projection, but geodetic distances calculated on a spheroid (Vincenty's formulae) or a sphere (Great-cicle distance). For instance, geopy is able to calculate them.

Here’s an example usage of Vincenty distance:

>>> from geopy.distance import vincenty
>>> newport_ri = (41.49008, -71.312796)
>>> cleveland_oh = (41.499498, -81.695391)
>>> print(vincenty(newport_ri, cleveland_oh).miles)
538.3904451566326

Using great-circle distance:

>>> from geopy.distance import great_circle
>>> newport_ri = (41.49008, -71.312796)
>>> cleveland_oh = (41.499498, -81.695391)
>>> print(great_circle(newport_ri, cleveland_oh).miles)
537.1485284062816

Source: https://pypi.python.org/pypi/geopy#measuring-distance


how to find point B with the above-described method?

You need to solve the direct problem, i.e. given an initial point, its azimuth and a geodesic distance calculate the final point. A Python implementation of the direct problem is available in the PyGeodesy package.

Example:

>>> from pygeodesy.ellipsoidalVincenty import LatLon 
>>> p = LatLon(-37.95103, 144.42487)
>>> d = p.destination(54972.271, 306.86816)
>>> print d.lon, d.lat
143.926497668 -37.6528177174
| improve this answer | |
  • Thanks @afalciano! This opens up a new perspective. Since I don't want to measure distance, I am wondering if you have a recommendation on how to find point B with the above-described method? – ustroetz Apr 23 '17 at 9:50
  • 1
    I suggest using the pygeodesy package to solve the direct problem. See the updated answer above. – Antonio Falciano Apr 23 '17 at 11:08
3

For this, I have to project point A to a projection that preserves distance (Equidistant).

You may have misunderstood "Equidistant" in Equidistant projections. It's a common misunderstanding that "Equidistant" projections preserves distances everywhere, which is impossible due to one of Gauss's theorems.

Some projections are called "Equidistant" because they preserve distance along a set of specific lines, such as meridians or parallels, which are called standard lines.

If in your case, points A and B are always to the East/West of each other, then in Sinusoidal projection, all parallel lines are standard.

Otherwise, you can use the spheroidal distance e.g. by computing distance directly with your points in lat/long (EPSG 4326) using PostGIS's geography type.

| improve this answer | |
  • Thank you @tinlyx for the detailed answer. I edited my question a bit by specifying the technology I am using (PostGIS is unfortunately not an option) and added the requirement that I want to create also a point to the north of point A (which I think eliminates the Sinusoidal projection as an option). – ustroetz Apr 23 '17 at 9:59
2

If you do not need accuracy down to the millimeter, you can just use one of the common transverse Mercator coordinate systems for this area, UTM zone 33 (EPSG 32633) or Gauß-Krüger zone 4 (EPSG 31468):

gk4_data = latlon_data.to_crs(epsg=31468)

Berlin is near a zone border in both systems, which results in small errors. As mentioned by Matte, in UTM, deviations are around 25 cm per km; Gauß-Krüger uses narrower zones, but without UTM's scale factor, so the deviations are approximately the same.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.