10

I'd like to create a (web)map, that is showing the range of an Airplane around some airports.

I tried to calculate a buffer with the range of the airplane. Here you can see the result here.

But now I have realized that the result is wrong because planes do not take the straight route but fly a curve because its shorter.

Is there a way to calculate the range with the shorter curve?

20

You could use the proj4 library to describe a circle using the great-circle distance.

For example, here's 3000km radius from Edinburgh, Tokyo, Cape Town and Quito in wgs84/Equirectangular. Only Quito is vaguely 'round', due to its proximity to the equator. I've also added in a single densified spoke line at an azimuth of 36 degrees (approx NE)

enter image description here

If we change to an azimuthal equidistant projection centred on Edinburgh, you'll see the radius around Edinburgh resolve to a circle...

enter image description here

On Mercator (like your web app), you see more distortion as you move away from the equator, but the buffers are more elliptical.

enter image description here

The following python code does that (requires pyproj and shapely)

import pyproj
from shapely.geometry import Polygon, MultiPoint, LineString
import math

def geodesicpointbuffer(longitude, latitude,
                        segments, distance_m,
                        geom_type=MultiPoint):
    """
    Creates a buffer in meters around a point given as long, lat in WGS84
    Uses the geodesic, so should be more accurate over larger distances

    :param longitude: center point longitude
    :param latitude: center point latitude
    :param segments: segments to approximate (more = smoother)
    :param distance_m: distance in meters
    :param geom_type: shapely type (e.g. Multipoint, Linestring, Polygon)
    :return: tuple (proj4 string, WKT of buffer geometry)
    """
    geodesic = pyproj.Geod(ellps='WGS84')
    coords = []
    for i in range(0, segments):
        angle = (360.0 / segments) * float(i)
        x1, y1, z1 = geodesic.fwd(lons=longitude,
                                  lats=latitude,
                                  az=angle,
                                  dist=distance_m,
                                  radians=False)
        coords.append((x1, y1))
        # makes a great circle for one spoke.
        if i==200:
            example = geodesic.npts(longitude,latitude,x1,y1,1000)
            coords2 = []
            for xx,yy in example:
                coords2.append((xx,yy))
            coords2.append((x1,y1)) # make sure we include endpoint ;-)
            flight = LineString(coords2)
            print(flight.wkt)

    ring = geom_type(coords)
    return "+init=EPSG:4326", ring.wkt


def main():
    # example : Cape Town. 3000km buffer.
    spec, wkt = geodesicpointbuffer(18.4637082653, -33.8496404007, 2000, 3000000.0, Polygon)
    print(spec)
    print(wkt)

if __name__ == "__main__":
    main()

You can paste the WKT output into QGIS using the useful QuickWKT plugin.

You could use other methods - as coneypylon mentioned, you could create a circle on a custom equidistant projection in meters, centred on your starting point. I find though that for large distances an error creeps in (only a few km at 2000 km, but for intercontinental distances these errors can mount up)

From memory, the mmqgis plugin allows buffering in km. I'm not sure which method it uses, though.

Note that you might have problems rendering polygons in QGIS that cross the antimeridian if you're starting in Asia - ogr2ogr with the -wrapdateline option can help here. You might find this is less of a problem with openlayers/leaflet, IIRC they allow longitudes greater than 180 and less than -180.

There's a good writeup about geodesic buffering here on the esri blog.

8

Depending on where your distance info is coming from, this may not matter. If you have a simple number giving the distance, the distance will be the same on any map projection which shows distance accurately (not Mercator, think pretty much any "equidistant" projection, such as an azimuthal orthographic projection or similar. A conformal projection, like Lambert Conformal Conic will do a reasonably okay job at the distance.). If you calculate and create the buffers in an equidistant projection, they will be (fairly) accurate, see here on how distance is calculated: ArcGIS Help

Be sure to set the layer's coordinate system is in an equidistant projection, not just the data frame.

Once calculated, the buffer will warp accordingly when put into Web Mercator or whatever other web projection you intend on using.

As far as why the lines themselves are curved, and why this might create problems:

The key problem is that plane routes on a Mercator projection such as this one are displayed as curved, like so:

Air Canada route map

This is a fundamental problem with Mercator maps, as they are intended for nautical navigation, where the properties of straight lines on these projections is valuable (a straight line on a Mercator projection is a rhumb line; a line with the same compass heading through the entire journey).

However, planes do not fly on rhumb lines, because fuel efficiency is more important than simple navigation, and thus fly along Great Circles, which appear as curves on a Mercator projection.

  • fuel efficiency and getting to the destination as quickly as possible. – cffk Apr 24 '17 at 21:16
  • ...which why the actual path will take the Jet Stream into account. – Vince Apr 25 '17 at 1:31
  • 2
    Gall-Peters is equal-area projection not equidistant. For equidistant you want something like an Azimuthal Orthographic Projection centered on your source. – HeikkiVesanto Apr 25 '17 at 8:35
  • Yeah, I was not thinking. A conformal projection will also do a reasonable job at preserving distance, yeah – coneypylon Apr 25 '17 at 14:12

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.