Map projections and their distortions

Draw India in Web MercatorWeb Mercator (EPSG:3857): The projection of web map tiles: the formulas of Mercator on a sphere applied to WGS 84 latitude and longitude. EPSG does not count it as a geodetic system, and it enlarges areas away from the equator., the projection of web maps, and a square kilometre of ground near Delhi measures 1.302 km² on the map while one near Chennai measures 1.060 km²: two equal fields, 23 % apart. In EPSG:7755, the CRS of Survey of India’s boundary data, the same two fields come out 3.2 % and 0.4 % small instead, with their shapes kept exactly.

No flat map can avoid such errors. This lesson shows which property each projection keeps, how to read the scale factors that say how much it stretches at a place, why UTM is reliable inside its zone and not beyond it, why squares drawn in UTM look turned on a web map, and how to choose a CRS for an India-wide study, a city’s roads and a web map.

Every flat map stretches something

Peel an orange and press the peel flat: it tears or stretches. A map does the same to the Earth. A map projectionmap projection: A set of formulas that turns latitude and longitude on an ellipsoid into flat x and y coordinates, such as eastings and northings in metres. Every projection stretches some lengths. turns latitude and longitude on an ellipsoid into flat x and y, and the QGIS manual states the consequence plainly: every map shows distortions of angular conformity, distance and area, and “it is usually impossible to preserve all characteristics at the same time in a map projection” (QGIS).

So each projection chooses what to keep. A conformal projectionconformal projection: A projection that keeps angles, and so the shapes of small features, at every point, such as Mercator, Lambert conic conformal and transverse Mercator. It cannot keep areas everywhere. keeps angles, and with them the shapes of small features. An equal-area projectionequal-area projection: A projection that keeps the areas of all features in proportion to their areas on the Earth, such as the Albers and Lambert azimuthal equal-area projections. It distorts angles and shapes. keeps areas in proportion. None keeps both, and none keeps every distance:

PROJ, the projection library that pyproj and PostGIS use, classifies the Lambert conic as “conformal conic” and notes its use in “many national and regional mapping systems”; it describes transverse Mercator as “a conformal projection in which a chosen meridian projects to a straight line at constant scale” (PROJ: Lambert conic, PROJ: transverse Mercator). The method’s name in a CRS definition therefore tells you what the CRS keeps.

Scale factors say how much, at each place

How much a projection stretches changes from place to place. A scale factorscale factor: A short length on a map grid divided by the same length on the ellipsoid, at one place and in one direction. Above 1 the grid makes lengths longer; below 1, shorter. measures it at one place and in one direction: a short length on the grid divided by the same length on the ellipsoid. The meridian scale h is the factor north-south, the parallel scale k east-west, and the areal scaleareal scale: A small area on a map grid divided by the same area on the ellipsoid, at one place. It is 1 in an equal-area projection and about 1.118 in Web Mercator at Pune. the factor for a small area. In a conformal projection h equals k, so the areal scale is k²; in an equal-area projection the areal scale is 1 and h and k differ.

The module below implements four projections, each with the formulas of its EPSG method, and measures the factors by stepping a millionth of a degree north and east of a point. Its first lines run the worked examples published with those formulas, so you can see that the code reproduces them:

"""Four map projections on an ellipsoid, and their scale factors, for the other programs of this lesson.

Formulas: IOGP report 373-07-2 (EPSG methods 1024, 9802 and 9807) and USGS Professional Paper 1395 (Albers).
"""
import math


class Ellipsoid:
    def __init__(self, a, inv_f):
        f = 1 / inv_f
        self.a, self.e2 = a, f * (2 - f)
        self.e, self.n = math.sqrt(self.e2), f / (2 - f)


WGS84 = Ellipsoid(6_378_137.0, 298.257223563)


class WebMercator:
    """EPSG:3857: the spherical Mercator formulas applied to WGS 84 latitude and longitude, with R = a."""

    def __init__(self, ell=WGS84):
        self.ell = ell

    def forward(self, lat, lon):
        a = self.ell.a
        return a * math.radians(lon), a * math.log(math.tan(math.pi / 4 + math.radians(lat) / 2))


class LambertConic:
    """Lambert Conic Conformal with two standard parallels (EPSG method 9802)."""

    def __init__(self, lat0, lon0, lat1, lat2, x0, y0, ell=WGS84):
        self.ell, self.lon0, self.x0, self.y0 = ell, lon0, x0, y0
        m1, m2 = self._m(lat1), self._m(lat2)
        t1, t2 = self._t(lat1), self._t(lat2)
        self.n = (math.log(m1) - math.log(m2)) / (math.log(t1) - math.log(t2))
        self.big_f = m1 / (self.n * t1 ** self.n)
        self.r0 = ell.a * self.big_f * self._t(lat0) ** self.n

    def _m(self, lat):
        p = math.radians(lat)
        return math.cos(p) / math.sqrt(1 - self.ell.e2 * math.sin(p) ** 2)

    def _t(self, lat):
        p, e = math.radians(lat), self.ell.e
        return math.tan(math.pi / 4 - p / 2) / ((1 - e * math.sin(p)) / (1 + e * math.sin(p))) ** (e / 2)

    def forward(self, lat, lon):
        r = self.ell.a * self.big_f * self._t(lat) ** self.n
        theta = self.n * math.radians(lon - self.lon0)
        return self.x0 + r * math.sin(theta), self.y0 + self.r0 - r * math.cos(theta)


class TransverseMercator:
    """Transverse Mercator (EPSG method 9807), with the 4th-order series of IOGP report 373-07-2."""

    def __init__(self, lat0, lon0, k0, x0, y0, ell=WGS84):
        self.ell, self.lon0, self.k0, self.x0, self.y0 = ell, lon0, k0, x0, y0
        n = ell.n
        self.b = ell.a / (1 + n) * (1 + n ** 2 / 4 + n ** 4 / 64)
        self.h = (n / 2 - 2 * n ** 2 / 3 + 5 * n ** 3 / 16 + 41 * n ** 4 / 180,
                  13 * n ** 2 / 48 - 3 * n ** 3 / 5 + 557 * n ** 4 / 1440,
                  61 * n ** 3 / 240 - 103 * n ** 4 / 140,
                  49561 * n ** 4 / 161280)
        self.m0 = 0.0 if lat0 == 0 else self.b * self._xi(self._beta(lat0), 0.0)[0]

    def _beta(self, lat):
        p, e = math.radians(lat), self.ell.e
        return math.atan(math.sinh(math.asinh(math.tan(p)) - e * math.atanh(e * math.sin(p))))

    def _xi(self, beta, dlon):
        eta0 = math.atanh(math.cos(beta) * math.sin(dlon))
        xi0 = math.asin(math.sin(beta) * math.cosh(eta0))
        xi = xi0 + sum(h * math.sin(2 * j * xi0) * math.cosh(2 * j * eta0) for j, h in enumerate(self.h, 1))
        eta = eta0 + sum(h * math.cos(2 * j * xi0) * math.sinh(2 * j * eta0) for j, h in enumerate(self.h, 1))
        return xi, eta

    def forward(self, lat, lon):
        xi, eta = self._xi(self._beta(lat), math.radians(lon - self.lon0))
        return self.x0 + self.k0 * self.b * eta, self.y0 + self.k0 * (self.b * xi - self.m0)

    def inverse(self, x, y):
        """Easting and northing back to latitude and longitude (degrees)."""
        n, e = self.ell.n, self.ell.e
        hi = (n / 2 - 2 * n ** 2 / 3 + 37 * n ** 3 / 96 - n ** 4 / 360,
              n ** 2 / 48 + n ** 3 / 15 - 437 * n ** 4 / 1440,
              17 * n ** 3 / 480 - 37 * n ** 4 / 840,
              4397 * n ** 4 / 161280)
        eta_ = (x - self.x0) / (self.b * self.k0)
        xi_ = ((y - self.y0) + self.k0 * self.m0) / (self.b * self.k0)
        xi0 = xi_ - sum(h * math.sin(2 * j * xi_) * math.cosh(2 * j * eta_) for j, h in enumerate(hi, 1))
        eta0 = eta_ - sum(h * math.cos(2 * j * xi_) * math.sinh(2 * j * eta_) for j, h in enumerate(hi, 1))
        beta = math.asin(math.sin(xi0) / math.cosh(eta0))
        q = q0 = math.asinh(math.tan(beta))
        for _ in range(15):
            q = q0 + e * math.atanh(e * math.tanh(q))
        lon = self.lon0 + math.degrees(math.asin(math.tanh(eta0) / math.cos(beta)))
        return math.degrees(math.atan(math.sinh(q))), lon


class AlbersEqualArea:
    """Albers Equal Area conic (EPSG method 9822), ellipsoidal formulas of USGS Professional Paper 1395."""

    def __init__(self, lat0, lon0, lat1, lat2, x0=0.0, y0=0.0, ell=WGS84):
        self.ell, self.lon0, self.x0, self.y0 = ell, lon0, x0, y0
        m1, m2 = self._m(lat1), self._m(lat2)
        q1, q2 = self._q(lat1), self._q(lat2)
        self.n = (m1 * m1 - m2 * m2) / (q2 - q1)
        self.c = m1 * m1 + self.n * q1
        self.rho0 = ell.a * math.sqrt(self.c - self.n * self._q(lat0)) / self.n

    def _m(self, lat):
        p = math.radians(lat)
        return math.cos(p) / math.sqrt(1 - self.ell.e2 * math.sin(p) ** 2)

    def _q(self, lat):
        s, e, e2 = math.sin(math.radians(lat)), self.ell.e, self.ell.e2
        return (1 - e2) * (s / (1 - e2 * s * s) - math.log((1 - e * s) / (1 + e * s)) / (2 * e))

    def forward(self, lat, lon):
        rho = self.ell.a * math.sqrt(self.c - self.n * self._q(lat)) / self.n
        theta = self.n * math.radians(lon - self.lon0)
        return self.x0 + rho * math.sin(theta), self.y0 + self.rho0 - rho * math.cos(theta)


def utm(zone):
    """WGS 84 / UTM zone (north): EPSG:326zz, central meridian 6 x zone - 183 degrees, scale 0.9996."""
    return TransverseMercator(0.0, 6 * zone - 183, 0.9996, 500_000.0, 0.0)


CRS = {
    "EPSG:3857 Web Mercator": WebMercator(),
    "EPSG:7755 India NSF LCC": LambertConic(24.0, 80.0, 12.472955, 35.1728044444444, 4_000_000.0, 4_000_000.0),
    "EPSG:32643 UTM 43N": utm(43),
    "EPSG:32644 UTM 44N": utm(44),
    "Albers equal-area for India": AlbersEqualArea(24.0, 80.0, 12.472955, 35.1728044444444),
}


def factors(proj, lat, lon, d=1e-6):
    """Scale factors at a place, by small steps north and east measured against the WGS 84 ellipsoid.

    Returns (h, k, areal scale, angular distortion in degrees, grid convergence in degrees): h along the meridian,
    k along the parallel; the convergence is the angle from true north to grid north, positive clockwise.
    """
    e2, a = proj.ell.e2, proj.ell.a
    p = math.radians(lat)
    w = math.sqrt(1 - e2 * math.sin(p) ** 2)
    rho, nu = a * (1 - e2) / w ** 3, a / w  # radii of curvature: along the meridian, along the prime vertical
    (xn, yn), (xs, ys) = proj.forward(lat + d, lon), proj.forward(lat - d, lon)
    (xe, ye), (xw, yw) = proj.forward(lat, lon + d), proj.forward(lat, lon - d)
    step = math.radians(2 * d)
    xp, yp, xl, yl = (xn - xs) / step, (yn - ys) / step, (xe - xw) / step, (ye - yw) / step
    h = math.hypot(xp, yp) / rho
    k = math.hypot(xl, yl) / (nu * math.cos(p))
    areal = (xl * yp - xp * yl) / (rho * nu * math.cos(p))
    big = math.sqrt(h * h + k * k + 2 * areal) / 2 + math.sqrt(max(h * h + k * k - 2 * areal, 0)) / 2
    small = math.sqrt(h * h + k * k + 2 * areal) / 2 - math.sqrt(max(h * h + k * k - 2 * areal, 0)) / 2
    omega = math.degrees(2 * math.asin((big - small) / (big + small)))
    return h, k, areal, omega, -math.degrees(math.atan2(xp, yp))


if __name__ == "__main__":
    # The worked examples of the formulas' sources, which these projections must reproduce.
    clarke = Ellipsoid(6_378_206.4, 294.9786982)
    us_ft = 1200 / 3937
    texas = LambertConic(27 + 50 / 60, -99.0, 28 + 23 / 60, 30 + 17 / 60, 2_000_000 * us_ft, 0.0, clarke)
    x, y = texas.forward(28.5, -96.0)
    print(f"Lambert conic, IOGP example (NAD27 / Texas South Central): {x / us_ft:,.2f} {y / us_ft:,.2f} US ft")
    x, y = utm(43).forward(18.5204, 73.8567)
    print(f"Transverse Mercator, Pune in UTM 43N (as GeographicLib's GeoConvert gives it): {x:,.2f} {y:,.2f} m")
    x, y = WebMercator().forward(24 + 22 / 60 + 54.433 / 3600, -(100 + 20 / 60))
    print(f"Web Mercator, IOGP example: {x:,.2f} {y:,.2f} m")
    x, y = AlbersEqualArea(23.0, -96.0, 29.5, 45.5, ell=clarke).forward(35.0, -75.0)
    print(f"Albers, USGS example (Clarke 1866): {x:,.1f} {y:,.1f} m")

Lambert conic, IOGP example (NAD27 / Texas South Central): 2,963,503.91 254,759.80 US ft
Transverse Mercator, Pune in UTM 43N (as GeographicLib's GeoConvert gives it): 379,320.82 2,048,144.62 m
Web Mercator, IOGP example: -11,169,055.58 2,800,000.00 m
Albers, USGS example (Clarke 1866): 1,885,472.7 1,535,925.0 m

Every value matches its source: the IOGP formulas’ worked examples to the hundredth of a foot or metre, the UTM coordinates of Pune as GeographicLib’s GeoConvert gives them (GeoConvert), and the USGS manual’s Albers example. The next program measures five CRSs at five cities:

"""Scale factors of five CRSs at five Indian cities, measured against the WGS 84 ellipsoid."""
from projections import CRS, factors

PLACES = {"Chennai": (13.0827, 80.2707), "Pune": (18.5204, 73.8567), "Kolkata": (22.5726, 88.3639),
          "Guwahati": (26.1445, 91.7362), "Delhi": (28.6139, 77.2090)}

for name, proj in CRS.items():
    print(f"{name}")
    print("  place       meridian  parallel  areal   angles  convergence")
    for place, (lat, lon) in PLACES.items():
        h, k, areal, omega, gamma = factors(proj, lat, lon)
        print(f"  {place:<10} {h:8.4f} {k:9.4f} {areal:7.4f} {omega:6.2f}° {gamma:9.2f}°")
print("1 km² of ground near Pune measures, in each CRS:")
for name, proj in CRS.items():
    print(f"  {name:<28} {factors(proj, *PLACES['Pune'])[2]:.3f} km²")

EPSG:3857 Web Mercator
  place       meridian  parallel  areal   angles  convergence
  Chennai      1.0330    1.0265  1.0604   0.37°     -0.00°
  Pune         1.0607    1.0543  1.1182   0.35°     -0.00°
  Kolkata      1.0886    1.0824  1.1784   0.33°     -0.00°
  Guwahati     1.1193    1.1133  1.2461   0.31°     -0.00°
  Delhi        1.1442    1.1383  1.3023   0.30°     -0.00°
EPSG:7755 India NSF LCC
  place       meridian  parallel  areal   angles  convergence
  Chennai      0.9980    0.9980  0.9960   0.00°      0.11°
  Pune         0.9849    0.9849  0.9701   0.00°     -2.50°
  Kolkata      0.9808    0.9808  0.9620   0.00°      3.40°
  Guwahati     0.9812    0.9812  0.9628   0.00°      4.77°
  Delhi        0.9838    0.9838  0.9678   0.00°     -1.13°
EPSG:32643 UTM 43N
  place       meridian  parallel  areal   angles  convergence
  Chennai      1.0037    1.0037  1.0073   0.00°      1.20°
  Pune         0.9998    0.9998  0.9996   0.00°     -0.36°
  Kolkata      1.0233    1.0233  1.0472   0.00°      5.21°
  Guwahati     1.0350    1.0350  1.0712   0.00°      7.55°
  Delhi        1.0002    1.0002  1.0004   0.00°      1.06°
EPSG:32644 UTM 44N
  place       meridian  parallel  areal   angles  convergence
  Chennai      0.9997    0.9997  0.9994   0.00°     -0.17°
  Pune         1.0067    1.0067  1.0134   0.00°     -2.28°
  Kolkata      1.0067    1.0067  1.0135   0.00°      2.84°
  Guwahati     1.0140    1.0140  1.0281   0.00°      4.78°
  Delhi        1.0013    1.0013  1.0026   0.00°     -1.82°
Albers equal-area for India
  place       meridian  parallel  areal   angles  convergence
  Chennai      1.0019    0.9981  1.0000   0.22°      0.11°
  Pune         1.0149    0.9854  1.0000   1.69°     -2.43°
  Kolkata      1.0194    0.9810  1.0000   2.20°      3.31°
  Guwahati     1.0193    0.9810  1.0000   2.19°      4.65°
  Delhi        1.0169    0.9834  1.0000   1.92°     -1.11°
1 km² of ground near Pune measures, in each CRS:
  EPSG:3857 Web Mercator       1.118 km²
  EPSG:7755 India NSF LCC      0.970 km²
  EPSG:32643 UTM 43N           1.000 km²
  EPSG:32644 UTM 44N           1.013 km²
  Albers equal-area for India  1.000 km²

The figure draws the areal scale column, the share by which each CRS enlarges or shrinks a small area:

The three conformal CRSs show h equal to k and no angular distortion, and still change areas: EPSG:7755 shrinks them by 0.4 % at Chennai and 3.8 % at Kolkata. The Albers conic, built here with the same standard parallels as EPSG:7755, keeps every area at 1.000 and pays with up to 2.2° of angular distortion. Web Mercator enlarges areas everywhere, by 6.0 % at Chennai and 30.2 % at Delhi.

Web Mercator is for display, not for measuring

Web Mercator, EPSG:3857, is the projection of web map tiles. EPSG describes it as using a “spherical development of ellipsoidal coordinates”: it applies the formulas of Mercator on a sphere to WGS 84 latitudes and longitudes, and EPSG lists it as “not a recognised geodetic system”, for web mapping and visualisation (EPSG:3857). PROJ notes that the computation uses the ellipsoid’s semi-major axis on a sphere, which makes the projection “non conformant if used with ellipsoid” (PROJ).

The factors program shows both effects. Web Mercator’s h and k differ slightly, by 0.30 to 0.37° of angular distortion, and its areal scale grows with latitude:

Areas measured in it are wrong by 6 to 30 % across these cities, and the error is not the same from district to district, so even a ranking by area comes out wrong.

UTM is accurate inside its zone and drifts outside it

UTMUTM (Universal Transverse Mercator): A system of 60 zones, each 6° of longitude wide and each a transverse Mercator projection with a scale factor of 0.9996 on its central meridian. Mainland India spans zones 42N to 47N. divides the world into 60 zones, each 6° of longitude wide, and gives each its own transverse Mercator projection. EPSG:32643, UTM zone 43N, is centred on the central meridiancentral meridian: The meridian on which a transverse Mercator or conic projection is centred, such as 75° E for UTM zone 43N. Along it, grid north is true north. 75° E with a scale factor of 0.9996 there, a false easting of 500,000 m, and an area of use from 72° E to 78° E (EPSG:32643). Mainland India, from about 68° E to 97° E, spans zones 42N to 47N, EPSG:32642 to EPSG:32647.

The scale factor of 0.9996 on the central meridian is deliberate: lengths there are 0.04 % short, the scale reaches 1 about 180 km either side, and it rises to about 1.0008 at the zone’s edge near 20° N. Inside its zone UTM is the most accurate grid in this lesson; outside it the error grows fast. The factors program measures zone 43N at the five cities:

Pune and Delhi lie inside the zone, and their scale is within 0.02 % of 1. Chennai, Kolkata and Guwahati lie in other zones, and zone 43N’s error there grows to 7.1 % in area; in its own zone, 44N, Chennai’s areal scale is 0.9994.

Grid north is not true north

On a projected grid, north is the direction in which the northing grows. True north is the direction along the meridian. The angle between the two, the grid convergencegrid convergence: The angle between true north, along the meridian, and grid north, the direction in which a projected CRS's northing grows, at one place. It is zero on the projection's central meridian., is zero on a projection’s central meridian and grows away from it. The factors program prints it in its last column; the figure shows three cases:

At Pune, 1.1° west of zone 43N’s central meridian, the grid’s north is 0.36° west of true north. At Guwahati, 11.7° east of EPSG:7755’s central meridian of 80° E, it is 4.77° east of true north, and in zone 43N, far outside the zone, 7.55°. Web Mercator draws meridians as vertical lines, so its grid north is true north everywhere. A feature aligned with a UTM grid therefore appears turned on a web map. The program builds a 2 km square ward in UTM zone 43N near Pune and draws it in EPSG:3857:

"""A square ward drawn in UTM zone 43N near Pune, as a north-up web map in EPSG:3857 shows it."""
import math

from projections import WebMercator, utm

zone, web = utm(43), WebMercator()
corners = [(378_000, 2_046_000), (380_000, 2_046_000), (380_000, 2_048_000), (378_000, 2_048_000)]  # 2 km x 2 km
lat_lon = [zone.inverse(x, y) for x, y in corners]
on_web = [web.forward(lat, lon) for lat, lon in lat_lon]
(x1, y1), (x2, y2) = on_web[0], on_web[1]
turn = math.degrees(math.atan2(y2 - y1, x2 - x1))
print(f"South-west corner: {lat_lon[0][0]:.5f}° N, {lat_lon[0][1]:.5f}° E")
print(f"The ward's south edge, due east in the UTM grid, points {turn:.3f}° anticlockwise of east on the web map")
side = math.hypot(x2 - x1, y2 - y1)
print(f"and is {side:,.0f} web-map metres long for 2,000 m of grid")

South-west corner: 18.50094° N, 73.84432° E
The ward's south edge, due east in the UTM grid, points 0.366° anticlockwise of east on the web map
and is 2,109 web-map metres long for 2,000 m of grid

The ward’s south edge, due east on the UTM grid, points 0.366° anticlockwise of east on the web map, and its 2,000 m are drawn 2,109 m long. Neither is an error in the data: they are the two grids’ convergence and scale, measured on the same square.

Choose the CRS from what you will measure

Every CRS above suits some jobs and spoils others, so start from the measurement. Answer the two questions in the table below about your own task; it rules out the CRSs that would distort what you measure and says why, with this lesson’s numbers:

In short: area statistics across India need an equal-area projection such as the Albers conic, or areas computed on the ellipsoid; lengths within one city need the city’s UTM zone or the ellipsoid; a web map draws in EPSG:3857 and measures nothing in it. EPSG:7755, conformal with two standard parallelsstandard parallel: A parallel of latitude along which a conic projection keeps lengths true, with a scale factor of exactly 1. EPSG:7755 has two, at about 12.47° N and 35.17° N. at about 12.47° N and 35.17° N, keeps shapes across the whole country and is a good common grid for national layers, as long as its areas are corrected or measured another way (EPSG:7755).

Practise: find the UTM zone of a point

Write the function that every UTM workflow starts with: given a point, the EPSG code of its zone. The sample tests use the lesson’s five cities, the edges of zones and the southern hemisphere:

def utm_epsg(lat, lon):
    """The EPSG code of the WGS 84 UTM zone for a point: 326zz in the north, 327zz in the south."""
    # Replace this line: work out the zone from the longitude, then add 32600 or 32700.
    return 32643

from zones import utm_epsg


def test_indian_cities():
    """finds the zone of five cities, from 43N in the west to 46N in the east"""
    assert utm_epsg(18.5204, 73.8567) == 32643  # Pune
    assert utm_epsg(28.6139, 77.2090) == 32643  # Delhi
    assert utm_epsg(13.0827, 80.2707) == 32644  # Chennai
    assert utm_epsg(22.5726, 88.3639) == 32645  # Kolkata
    assert utm_epsg(26.1445, 91.7362) == 32646  # Guwahati


def test_zone_edges():
    """puts a point on a zone's western edge in that zone, and 180° in zone 60"""
    assert utm_epsg(20.0, 72.0) == 32643
    assert utm_epsg(20.0, 78.0) == 32644
    assert utm_epsg(10.0, 180.0) == 32660
    assert utm_epsg(10.0, -180.0) == 32601


def test_southern_hemisphere():
    """uses 327zz south of the equator and 326zz on it"""
    assert utm_epsg(-22.9068, -43.1729) == 32723  # Rio de Janeiro
    assert utm_epsg(0.0, 73.5) == 32643


def test_outside_utm():
    """refuses latitudes beyond 84° N and 80° S, where UTM does not apply"""
    for lat in (85.0, -81.0):
        try:
            utm_epsg(lat, 10.0)
        except ValueError:
            continue
        assert False, f"latitude {lat} should raise ValueError"

Then practise the six ideas, from what every flat map stretches to the choice of a CRS, beginning with the one your answers so far leave weakest:

Interview questions

Warm-up (fresher): Why can’t a flat map keep both areas and shapes?

Because the ellipsoid is curved in two directions at once and a plane is not, so flattening it must stretch some lengths. A projection can choose how: a conformal projection stretches equally in every direction at each point, which keeps angles and the shapes of small features, but the amount of stretch changes from place to place, so areas change; EPSG:7755 shrinks areas by 0.4 % at Chennai and 3.8 % at Kolkata. An equal-area projection balances a stretch in one direction with a shrink in the other, which keeps areas but bends angles; the Albers conic of this lesson keeps every area and distorts angles by up to 2.2°. Keeping both would need no stretching at all, which no flat map can do.

Two questions in the module’s bank take this further, one about ranking districts by area and one about a city that grows past its UTM zone. Each lists what it tests and a first step:

Key takeaways

• Every flat map stretches the Earth: a conformal projection keeps angles and small shapes, an equal-area projection keeps areas, and none keeps both.
• Scale factors measure the stretch at a place, h along the meridian, k along the parallel and the areal scale for areas; for Web Mercator the areal scale runs from 1.060 at Chennai to 1.302 at Delhi.
• Web Mercator is a spherical development for web tiles, not a geodetic system, so areas measured in it are 6 to 30 % too large across India.
• UTM is accurate inside its 6° zone, with a scale of 0.9998 at Pune in zone 43N, and drifts outside it, to an areal scale of 1.071 at Guwahati.
• Grid convergence, the angle between grid north and true north, is 0.36° at Pune in UTM 43N, so UTM squares look turned on a north-up web map.
• Choose the CRS from the measurement: equal-area or the ellipsoid for national areas, the local UTM zone for city lengths, and EPSG:3857 only for drawing.