Map Range of Airplanes

Question

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?

Answer 1

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)

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

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

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.

Answer 2

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:

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.