1

For each coordinate I have, I find the distance from the equator in kilometers giving me two distances: 

from pyproj import Geod
wgs84_geod = Geod(ellps='WGS84')    
_,_, lon_dist = wgs84_geod.inv(0, 0,lon, 0)
_,_, lat_dist = wgs84_geod.inv(0, 0,0, lat)

As a sanity check,I can recalculate the original coordinate from these values as follows (assume the direction from the equator coordinate (0,0) is North and West:

_, new_lat, _ = wgs84_geod.fwd(0,0, 0, lat_dist)
new_lon, _, _ = wgs84_geod.fwd(0, 0, 90, lon_dist)

This gives me back the same coordinates I started with.

Now I want to find the closest kilometer point to my coordinate. I round the lon_dist and lat_dist to kilometers from the equator values.

lat_km_dist = round(lat_dist/1000)*1000 #to nearest km and back to meters
lon_km_dist = round(lon_dist/1000)*1000

I get coordinates using these distances in the same way as before

_, km_lat, _ = wgs84_geod.fwd(0,0, 0, lat_km_dist)
km_lon, _, _ = wgs84_geod.fwd(0, 0, 90, lon_km_dist)

The logic should be that for multiple coordinates in the same area, the closest distance between any km_lat, km_lon pair should be 1km. This is true in the North/South axis, but for longitudes the distance varies depending on which latitude I'm at. I'm attaching two screenshots to visualize the problem where the km_lat, km_lon coordinates are represented by black circles at the center of polygons with an area of 1km.

How can I correct for this?

Copenhagen

San Francisco

clurhur
  • 81
  • 6

1 Answers1

1

What this algorithm is essentially doing is that it constructs an equidistant mesh (with points 1km apart) on the equator (lat=0) and the main meridian (lon=0). It then effectively constructs a grid on the ellipsoid as a Cartesian product of these points.

However, the lat/lon coordinates do not form a Cartesian frame, the resulting parallels/meridians generated by these grid points define "squares" the size of which depends not only on the particular longitude, but latitude as well. On a perfect sphere, this would work in the north-south direction since then an equidistant (in terms of great-circle distance) grid on lon=0 is also equidistant in the latitude (difference in latitude being equal to the difference in distance over the radius of the sphere).

In other words, if you fix two latitudes lat1, lat2 and for a particular longitude lon move from (lat1, lon), (lat2, lon) 1km in, say, westward direction, then these newly obtained points won't have the same longitude...

I am not fully sure what you are trying to achieve but if the goal is to obtain some representative points not too close to each other, then perhaps hierarchical clustering in terms of the great-circle distance could provide reasonable results...

EDIT:

As an approximative workaround, you could get most likely away by choosing another reference point than (0, 0) - the new reference point shouldn't be too far away from the area that you are trying to describe (something like a "bottom-left" corner of the area of interest). If the entire area of interest doesn't cover a significant part of the globe (large span of latitudes) then the discrepancies will be quite small so that they will be probably almost invisible in the GoogleMaps visualization...

So if you are interested in Denmark (judging by the screenshots), then something like the following might work:

lat_ref, lon_ref = 53.637976, 6.694138

_,_, lon_dist = wgs84_geod.inv(lon_ref,lat_ref, lon, 0)
_,_, lat_dist = wgs84_geod.inv(lon_ref,lat_ref, 0, lat)

lat_km_dist = round(lat_dist/1000)*1000 #to nearest km and back to meters
lon_km_dist = round(lon_dist/1000)*1000


_, km_lat, _ = wgs84_geod.fwd(lon_ref,lat_ref,  0, lat_km_dist)
km_lon, _, _ = wgs84_geod.fwd(lon_ref,lat_ref, 90, lon_km_dist)
ewcz
  • 12,819
  • 1
  • 25
  • 47
  • The goal is, as you say, to reduce a large number of coordinates to fewer, representative points. There is a need for this to be independent of the dataset, so given different coordinates in the same area, the underlying grid would not change. For this reason, I don't want to use clustering. Is there a way to factor in the effect of latitude on longitude distance in this method? – clurhur Mar 08 '17 at 14:32
  • @clurhur you might try to change the reference point - I have included an example... – ewcz Mar 08 '17 at 14:33
  • thanks for this. I have tried this approach and it gets me closer to the 1km distance I want. The difficulty is this relies on knowing these reference coordinates. As I want to use this method anywhere in the world, this requires gathering many such references. I have seen a list of country centers, but for countries that span a large North/South area, the same problem occurs. – clurhur Mar 08 '17 at 14:40
  • @clurhur so perhaps something like this then? http://stackoverflow.com/questions/9600801/evenly-distributing-n-points-on-a-sphere – ewcz Mar 08 '17 at 14:53
  • Thanks for the tips @ewcz. I have gone with the approach of using one of the original coordinates as the center instead of the origin (0,0). – clurhur Mar 13 '17 at 08:05