GIS & Geospatial Data (QGIS, GeoPandas, PostGIS) Module 2 – Coordinate reference systems
Measuring distance and area correctly
Measure distances and areas on the ellipsoid, predict the error of degrees, Web Mercator, a conic or the wrong axis order, and convert hectares and acres.
What you will learn
- Compute geodesic distances and areas on the WGS 84 ellipsoid, and say when a planar measurement is close enough
- Predict the error of measuring in degrees, in Web Mercator, in a conformal conic or in the local UTM zone
- Use QGIS $area and area(), PostGIS geography and DuckDB's spheroid functions with the right inputs
- Convert areas into hectares, acres and local land units from their stated definitions
Before you start
What you’ll learn · Before you start
What you will learn
- Compute geodesic distances and areas on the WGS 84 ellipsoid, and say when a planar measurement is close enough
- Predict the error of measuring in degrees, in Web Mercator, in a conformal conic or in the local UTM zone
- Use QGIS $area and area(), PostGIS geography and DuckDB's spheroid functions with the right inputs
- Convert areas into hectares, acres and local land units from their stated definitions
Before you start
On this page
From Pune to Mumbai is 120,138.17 m on the WGS 84 ellipsoid. Measure the same line as a straight segment in Web Mercator and you get 127,065 m, 5.8 % too long; in EPSG:7755, 1.5 % too short; pass latitude and longitude the wrong way round to a geodesic function and it reports 110,690 m. Areas go the same way: a 400 ha ward near Pune measures 447 ha in Web Mercator.
This lesson measures distances and areas on the ellipsoid, shows when a flat grid or a sphere is close enough, puts numbers on the error of each wrong method, shows how QGIS, PostGIS and DuckDB decide what they measure, and converts areas into hectares, acres and local land units.
A distance on the Earth is a geodesic
The shortest line between two points on the ellipsoid is a geodesic: The shortest line between two points on the ellipsoid; its length is the distance between them on the Earth's model, 120,138.17 m from Pune to Mumbai on WGS 84., and its length is the distance a survey or a map needs. Karney’s algorithms, implemented in GeographicLib, compute it with an error “less than 15 nanometers”; the older formula of Vincenty is accurate to about 0.1 mm but “fails to converge for nearly antipodal points” (GeographicLib, Karney). The programs below use Vincenty’s formula, which is exact enough for any two places in India, and the first one measures Pune to Mumbai seven ways:
"""Pune to Mumbai, measured seven ways."""
from measure import R_MEAN, geodesic, haversine, india_lcc, utm43, web_mercator
pune, mumbai = (18.5204, 73.8567), (19.0760, 72.8777) # latitude, longitude on WGS 84
true = geodesic(*pune, *mumbai)[0]
print(f"On the WGS 84 ellipsoid (geodesic): {true:,.2f} m")
def show(label, metres):
print(f" {label:<38} {metres:12,.2f} m {metres - true:+9,.2f} m {(metres - true) / true:+8.3%}")
show(f"on a sphere, R = {R_MEAN:,.1f} m", haversine(*pune, *mumbai))
for label, project in [("straight line in UTM 43N", utm43), ("straight line in EPSG:7755", india_lcc),
("straight line in EPSG:3857", web_mercator)]:
(x1, y1), (x2, y2) = project(*pune), project(*mumbai)
show(label, ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** 0.5)
show("geodesic with the axes swapped", geodesic(pune[1], pune[0], mumbai[1], mumbai[0])[0])
degrees = ((mumbai[0] - pune[0]) ** 2 + (mumbai[1] - pune[1]) ** 2) ** 0.5
print(f" straight line in degrees: {degrees:.4f}, a number with no unit of length")Output
On the WGS 84 ellipsoid (geodesic): 120,138.17 m on a sphere, R = 6,371,008.8 m 120,152.47 m +14.30 m +0.012% straight line in UTM 43N 120,135.35 m -2.82 m -0.002% straight line in EPSG:7755 118,276.90 m -1,861.27 m -1.549% straight line in EPSG:3857 127,065.38 m +6,927.21 m +5.766% geodesic with the axes swapped 110,690.41 m -9,447.76 m -7.864% straight line in degrees: 1.1257, a number with no unit of length
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 distances.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
Its formulas come from a module shared by both programs of the lesson: Vincenty’s inverse and direct formulas, the haversine formula, the three projections of the comparison and an equal-area projection for areas:
"""Distances and areas on the WGS 84 ellipsoid, on a sphere and in four projections, for this lesson's programs."""
import math
A, F = 6_378_137.0, 1 / 298.257223563 # WGS 84
B = A * (1 - F)
E2 = F * (2 - F)
E = math.sqrt(E2)
R_MEAN = (2 * A + B) / 3 # the mean radius of the ellipsoid, 6,371,008.8 m
def geodesic(lat1, lon1, lat2, lon2):
"""Vincenty's inverse formula: metres between two points on WGS 84, and the azimuth at the first (degrees)."""
u1, u2 = math.atan((1 - F) * math.tan(math.radians(lat1))), math.atan((1 - F) * math.tan(math.radians(lat2)))
lam = big_l = math.radians(lon2 - lon1)
for _ in range(200):
sin_s = math.hypot(math.cos(u2) * math.sin(lam),
math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lam))
if sin_s == 0:
return 0.0, 0.0
cos_s = math.sin(u1) * math.sin(u2) + math.cos(u1) * math.cos(u2) * math.cos(lam)
sigma = math.atan2(sin_s, cos_s)
sin_al = math.cos(u1) * math.cos(u2) * math.sin(lam) / sin_s
cos2 = 1 - sin_al ** 2
c2m = cos_s - 2 * math.sin(u1) * math.sin(u2) / cos2 if cos2 else 0.0
c = F / 16 * cos2 * (4 + F * (4 - 3 * cos2))
last = lam
lam = big_l + (1 - c) * F * sin_al * (sigma + c * sin_s * (c2m + c * cos_s * (-1 + 2 * c2m ** 2)))
if abs(lam - last) < 1e-12:
break
u_sq = cos2 * (A * A - B * B) / (B * B)
big_a = 1 + u_sq / 16384 * (4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq)))
big_b = u_sq / 1024 * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))
d_s = big_b * sin_s * (c2m + big_b / 4 * (cos_s * (-1 + 2 * c2m ** 2)
- big_b / 6 * c2m * (-3 + 4 * sin_s ** 2) * (-3 + 4 * c2m ** 2)))
cross = math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lam)
az = math.atan2(math.cos(u2) * math.sin(lam), cross)
return B * big_a * (sigma - d_s), math.degrees(az)
def destination(lat, lon, azimuth, s):
"""Vincenty's direct formula: the point s metres from (lat, lon) along the geodesic that leaves at azimuth."""
al1 = math.radians(azimuth)
tan_u1 = (1 - F) * math.tan(math.radians(lat))
cos_u1 = 1 / math.sqrt(1 + tan_u1 ** 2)
sin_u1 = tan_u1 * cos_u1
sigma1 = math.atan2(tan_u1, math.cos(al1))
sin_al = cos_u1 * math.sin(al1)
cos2 = 1 - sin_al ** 2
u_sq = cos2 * (A * A - B * B) / (B * B)
big_a = 1 + u_sq / 16384 * (4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq)))
big_b = u_sq / 1024 * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))
sigma = s / (B * big_a)
for _ in range(200):
c2m = math.cos(2 * sigma1 + sigma)
d_s = big_b * math.sin(sigma) * (c2m + big_b / 4 * (math.cos(sigma) * (-1 + 2 * c2m ** 2) - big_b / 6 * c2m
* (-3 + 4 * math.sin(sigma) ** 2) * (-3 + 4 * c2m ** 2)))
last, sigma = sigma, s / (B * big_a) + d_s
if abs(sigma - last) < 1e-13:
break
c2m = math.cos(2 * sigma1 + sigma)
x = sin_u1 * math.sin(sigma) - cos_u1 * math.cos(sigma) * math.cos(al1)
up = sin_u1 * math.cos(sigma) + cos_u1 * math.sin(sigma) * math.cos(al1)
lat2 = math.atan2(up, (1 - F) * math.hypot(sin_al, x))
lam = math.atan2(math.sin(sigma) * math.sin(al1),
cos_u1 * math.cos(sigma) - sin_u1 * math.sin(sigma) * math.cos(al1))
c = F / 16 * cos2 * (4 + F * (4 - 3 * cos2))
big_l = lam - (1 - c) * F * sin_al * (sigma + c * math.sin(sigma)
* (c2m + c * math.cos(sigma) * (-1 + 2 * c2m ** 2)))
return math.degrees(lat2), lon + math.degrees(big_l)
def haversine(lat1, lon1, lat2, lon2, r=R_MEAN):
"""Great-circle metres on a sphere of radius r."""
p1, p2 = math.radians(lat1), math.radians(lat2)
h = math.sin((p2 - p1) / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(math.radians(lon2 - lon1) / 2) ** 2
return 2 * r * math.asin(math.sqrt(h))
def utm43(lat, lon):
"""WGS 84 / UTM zone 43N (EPSG:32643), with the 4th-order series of IOGP report 373-07-2."""
n = F / (2 - F)
b = A / (1 + n) * (1 + n ** 2 / 4 + n ** 4 / 64)
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)
p = math.radians(lat)
beta = math.atan(math.sinh(math.asinh(math.tan(p)) - E * math.atanh(E * math.sin(p))))
eta0 = math.atanh(math.cos(beta) * math.sin(math.radians(lon - 75)))
xi0 = math.asin(math.sin(beta) * math.cosh(eta0))
xi = xi0 + sum(c * math.sin(2 * j * xi0) * math.cosh(2 * j * eta0) for j, c in enumerate(h, 1))
eta = eta0 + sum(c * math.cos(2 * j * xi0) * math.sinh(2 * j * eta0) for j, c in enumerate(h, 1))
return 500_000 + 0.9996 * b * eta, 0.9996 * b * xi
def utm43_inverse(x, y):
"""Easting and northing in UTM zone 43N back to latitude and longitude on WGS 84."""
n = F / (2 - F)
b = A / (1 + n) * (1 + n ** 2 / 4 + n ** 4 / 64)
h = (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_, xi_ = (x - 500_000) / (0.9996 * b), y / (0.9996 * b)
xi0 = xi_ - sum(c * math.sin(2 * j * xi_) * math.cosh(2 * j * eta_) for j, c in enumerate(h, 1))
eta0 = eta_ - sum(c * math.cos(2 * j * xi_) * math.sinh(2 * j * eta_) for j, c in enumerate(h, 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))
return math.degrees(math.atan(math.sinh(q))), 75 + math.degrees(math.asin(math.tanh(eta0) / math.cos(beta)))
def _t(p):
return math.tan(math.pi / 4 - p / 2) / ((1 - E * math.sin(p)) / (1 + E * math.sin(p))) ** (E / 2)
def _m(p):
return math.cos(p) / math.sqrt(1 - E2 * math.sin(p) ** 2)
def india_lcc(lat, lon):
"""WGS 84 / India NSF LCC (EPSG:7755), Lambert Conic Conformal (2SP)."""
p0, p1, p2 = math.radians(24.0), math.radians(12.472955), math.radians(35.1728044444444)
n = (math.log(_m(p1)) - math.log(_m(p2))) / (math.log(_t(p1)) - math.log(_t(p2)))
f = _m(p1) / (n * _t(p1) ** n)
r, r0 = A * f * _t(math.radians(lat)) ** n, A * f * _t(p0) ** n
theta = n * math.radians(lon - 80.0)
return 4_000_000 + r * math.sin(theta), 4_000_000 + r0 - r * math.cos(theta)
def web_mercator(lat, lon):
"""WGS 84 / Pseudo-Mercator (EPSG:3857)."""
return A * math.radians(lon), A * math.log(math.tan(math.pi / 4 + math.radians(lat) / 2))
def _q(p):
s = math.sin(p)
return (1 - E2) * (s / (1 - E2 * s * s) - math.log((1 - E * s) / (1 + E * s)) / (2 * E))
def equal_area(lat, lon, lat0, lon0):
"""Lambert cylindrical equal-area, ellipsoidal, true scale at lat0: keeps every area of the ellipsoid."""
k0 = _m(math.radians(lat0))
return A * k0 * math.radians(lon - lon0), A * _q(math.radians(lat)) / (2 * k0)
def planar_area(points):
"""The shoelace formula: the area of a polygon from its vertices in a plane."""
return abs(sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(points, points[1:] + points[:1]))) / 2
def densify(vertices, step=1_000.0):
"""The polygon's edges as geodesics, with a vertex at least every `step` metres: a list of (lat, lon)."""
out = []
for (lat1, lon1), (lat2, lon2) in zip(vertices, vertices[1:] + vertices[:1]):
s, az = geodesic(lat1, lon1, lat2, lon2)
pieces = max(1, math.ceil(s / step))
out += [(lat1, lon1)] + [destination(lat1, lon1, az, s * i / pieces) for i in range(1, pieces)]
return out
def geodesic_area(vertices, step=2.0):
"""The area on WGS 84 of a polygon with geodesic edges, in m², measured in an equal-area projection."""
ring = densify(vertices, step)
lat0 = sum(lat for lat, _ in vertices) / len(vertices)
lon0 = sum(lon for _, lon in vertices) / len(vertices)
return planar_area([equal_area(lat, lon, lat0, lon0) for lat, lon in ring])The geodesic is 120,138.17 m; GeographicLib’s own online solver gives 120,138.170 m for the same two points
(GeodSolve). A sphere of the ellipsoid’s mean radius,
6,371,008.8 m, with the haversine formula: A formula for the great-circle distance between two points on a sphere from their latitudes and longitudes. With the Earth's mean radius it is fast and about 0.01 % off over distances such as Pune to Mumbai., gives 14.30 m more, 0.012 %. That is the
trade PostGIS offers too: ST_DistanceSphere uses a sphere of radius 6,371,008 m by default and is “faster than
ST_DistanceSpheroid, but less accurate” (PostGIS). For most
analysis within India a sphere is close enough; for surveys, legal distances or any check of others’ numbers, use the
geodesic.
Measure your own pair of places
Press Edit and change pune and mumbai to two places you know, latitude first. Compare the sphere’s error with
0.012 % for lines in other directions and at other latitudes.
A planar distance is only as good as the grid’s scale
The same program measures the line as a straight segment in three projected CRSs, the way a GIS measures geometry in the layer’s own grid. The chart puts every method against the geodesic:
Pune to Mumbai measured five ways, against the geodesic on WGS 84
Text description of the diagram
A chart of five horizontal bars around a vertical line that stands for the geodesic distance on the WGS 84 ellipsoid, 120,138.17 m, so an error of 0 %. A bar to the right means the method gives too long a distance; to the left, too short.
A sphere with the mean radius of 6,371,008.8 m (the haversine formula): 14.3 m too long, a bar too short to draw. A straight line in UTM zone 43N: 2.8 m too short, also too short to draw. A straight line in EPSG:7755: 1.55 % too short. A straight line in EPSG:3857, Web Mercator: 5.77 % too long, the longest bar to the right. The geodesic formula with latitude and longitude swapped: 7.86 % too short, the longest bar to the left.
In UTM zone 43N the straight line is 2.82 m short, 0.002 %; in EPSG:7755 it is 1.549 % short; in Web Mercator 5.766 % too long. Those are the scale factors of the previous lessons at work: near Pune, UTM 43N’s scale is 0.9998, EPSG:7755’s 0.9849 and Web Mercator’s about 1.054 to 1.061, and a planar length is wrong by the scale along the line, minus one. A straight line in degrees is 1.1257, which is no distance at all, because at India’s latitudes a degree of longitude is shorter than a degree of latitude, and it shrinks further towards the poles.
So a planar distance is close enough when the grid’s scale is close to 1 along the whole line: inside a city in its own UTM zone, the error here was under 3 m in 120 km. Across states, through several zones, or in any CRS made for display, measure the geodesic.
Key idea
A planar length carries the grid’s scale error: a few parts in a hundred thousand in the local UTM zone, 1.5 % in EPSG:7755 and 5.8 % in Web Mercator near Pune. Measure on the ellipsoid when the scale is not close to 1.
Checkpoint
Check: Geodesic distances · Planar distances and the grid's scale
Areas belong on the ellipsoid or in an equal-area grid
Areas follow the areal scale, which errs twice as much as the length scale in a conformal grid. The next program takes the ward of the projections lesson, a 2 km square in UTM zone 43N near Pune, and measures it on the ellipsoid and in three grids; then it does the same for a block of central India from 15° N to 25° N and 75° E to 85° E:
"""A 2 km square ward near Pune, and a 10° block of India, measured in several ways."""
from measure import densify, geodesic_area, india_lcc, planar_area, utm43, utm43_inverse, web_mercator
ward_utm = [(378_000, 2_046_000), (380_000, 2_046_000), (380_000, 2_048_000), (378_000, 2_048_000)]
ward = [utm43_inverse(x, y) for x, y in ward_utm] # the corners as latitude and longitude on WGS 84
true = geodesic_area(ward)
print(f"The ward, a 2 km x 2 km square in UTM zone 43N near Pune:")
print(f" on the ellipsoid (geodesic edges) {true:14,.2f} m²")
for label, area in [("planar, in UTM 43N", planar_area(ward_utm)),
("planar, in EPSG:7755", planar_area([india_lcc(*p) for p in ward])),
("planar, in EPSG:3857", planar_area([web_mercator(*p) for p in ward]))]:
print(f" {label:<35} {area:14,.2f} m² {(area - true) / true:+8.3%}")
print(f" that is {true / 10_000:.4f} ha or {true / 4_046.8564224:.4f} acres")
block = [(15.0, 75.0), (15.0, 85.0), (25.0, 85.0), (25.0, 75.0)] # latitude, longitude
true = geodesic_area(block, 500.0)
ring = densify(block, 500.0)
print(f"A block from 15° N to 25° N and 75° E to 85° E: {true / 1e6:,.0f} km² on the ellipsoid")
for label, project in [("EPSG:7755", india_lcc), ("EPSG:3857", web_mercator)]:
area = planar_area([project(*p) for p in ring])
print(f" planar in {label}: {area / 1e6:,.0f} km² {(area - true) / true:+.2%}")Output
The ward, a 2 km x 2 km square in UTM zone 43N near Pune: on the ellipsoid (geodesic edges) 4,001,752.60 m² planar, in UTM 43N 4,000,000.00 m² -0.044% planar, in EPSG:7755 3,882,247.90 m² -2.986% planar, in EPSG:3857 4,474,235.76 m² +11.807% that is 400.1753 ha or 988.8546 acres A block from 15° N to 25° N and 75° E to 85° E: 1,158,951 km² on the ellipsoid planar in EPSG:7755: 1,122,429 km² -3.15% planar in EPSG:3857: 1,323,828 km² +14.23%
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 areas.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
On the ellipsoid the ward has 4,001,752.60 m²; GeographicLib’s online planimeter gives 4,001,752.5 m² for its four corners (Planimeter). In its own UTM zone the planar area, exactly 4,000,000 m², is 0.044 % small; in EPSG:7755 it is 2.986 % small and in Web Mercator 11.807 % large. For the block, 1,158,951 km² on the ellipsoid, EPSG:7755 gives 3.15 % less and Web Mercator 14.23 % more. The program measures the ellipsoidal area by projecting the polygon, its edges divided into short geodesic pieces, onto an equal-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. More in “Map projections and their distortions”, which keeps every area of the ellipsoid.
Libraries warn about the worst case. GeoPandas computes areas “expressed in the units of the CRS”, notes that “area may be invalid for a geographic CRS using degrees as units”, and advises projecting first (GeoPandas). A projected CRS is not enough on its own: it must be equal-area, or local enough that its areal scale is close to 1.
Every engine measures what you ask it to
Each engine measures in two ways, and the name of the function or the type of the data decides which:
| Engine | On the ellipsoid | Planar, in the data’s own CRS |
|---|---|---|
| QGIS | $area and $length, with the project’s ellipsoid set |
area() and length(), “always planimetric” in the geometry’s SRS |
| PostGIS | geography columns: ST_Distance on the spheroid by default, a sphere with use_spheroid false |
geometry columns: ST_Distance in the units of the SRS |
| DuckDB spatial | ST_Distance_Spheroid, ST_Area_Spheroid: latitude, longitude order |
ST_Distance, the planar distance |
The table follows each engine’s documentation: QGIS for its expressions
(QGIS), PostGIS for ST_Distance,
which since version 2.2 uses GeographicLib on the spheroid (PostGIS), and
DuckDB, whose ST_Distance_Spheroid uses GeographicLib and whose spheroid functions expect “a [latitude, longitude]
axis order”
(DuckDB).
That last rule is a trap for data stored longitude first. Points built as ST_Point(longitude, latitude) and passed to
ST_Distance_Spheroid are read with their axes swapped, and the program shows what that does: 110,690.41 m instead of
120,138.17 m, 7.864 % short, for the swapped points GeographicLib’s solver also measures at 110,690.413 m. Answer the
three questions of the table below for your own measurement; it rules out the methods that would distort it:
Answer the questions above: the table then marks what fits your case and says why the rest does not.
| Method | Lengths in a city | Lengths across India | Areas | Millions of pairs |
|---|---|---|---|---|
| Geodesic on WGS 84 Works on: latitude, longitude | Good exact: 120,138.17 m from Pune to Mumbai | Good exact at any distance | Good exact: the ward is 4,001,752.6 m² | Mixed an iterative formula, slower than a sphere |
| Sphere (haversine) Works on: latitude, longitude | Good 14.3 m too long over 120 km, 0.012 % | Mixed small errors that depend on direction and latitude | Mixed spherical areas differ slightly from the ellipsoid | Good one closed formula per pair |
| Planar in the local UTM zone Works on: UTM metres | Good 2.8 m short over 120 km near Pune | Poor one zone spans 6°; outside it the scale grows | Good the ward 0.044 % small in zone 43N | Good plain arithmetic once projected |
| Planar in EPSG:7755 Works on: national LCC metres | Poor 1.55 % short from Pune to Mumbai | Mixed conformal, but 0.2 to 1.9 % short at the lesson's cities | Poor the ward 2.99 % small, a 10° block 3.15 % small | Good plain arithmetic once projected |
| Planar in EPSG:3857 Works on: web map metres | Poor 5.77 % too long from Pune to Mumbai | Poor errors grow with latitude | Poor the ward 11.81 % large, a 10° block 14.23 % large | Good plain arithmetic once projected |
| Planar in degrees Works on: raw latitude, longitude | Poor 1.1257, a number with no unit of length | Poor a degree of longitude shrinks with latitude | Poor square degrees, not square metres | Mixed fast and meaningless |
How the answers decide
- “What do you measure?” Areas: rules out Planar in EPSG:3857 and Planar in degrees and adds a caveat to Planar in EPSG:7755 and Sphere (haversine). Areas need an areal scale of 1: Web Mercator made the ward 11.81 % large, EPSG:7755 2.99 % small, degrees give square degrees, and a sphere's areas differ slightly from the ellipsoid's.
- “What do you measure?” Areas and “Where does it lie?” Across several zones or all India: rules out Planar in the local UTM zone and Planar in EPSG:7755. Across India no single conformal grid keeps areas; measure geodesic areas on the ellipsoid.
- “What do you measure?” Distances or lengths and “Where does it lie?” Across several zones or all India: rules out Planar in the local UTM zone, Planar in EPSG:3857 and Planar in degrees and adds a caveat to Planar in EPSG:7755 and Sphere (haversine). Long lines cross the scale changes of any flat grid; geodesics on the ellipsoid are exact at every length.
- “What do you measure?” Distances or lengths and “Where does it lie?” Inside one city or one UTM zone: rules out Planar in EPSG:3857 and Planar in degrees and adds a caveat to Planar in EPSG:7755. In one city the local UTM zone was 2.8 m short over 120 km, while EPSG:3857 was 5.77 % too long.
- “How many measurements at a time?” Millions of pairs, at speed: adds a caveat to Geodesic on WGS 84. For millions of pairs, a sphere or a local UTM grid is faster; check that their 0.012 % or zone-edge errors are acceptable.
Checkpoint
Check: Measuring areas · What each engine measures
Hectares and acres are exact; local units need their definition
Once an area is right in square metres, its units are a matter of exact definitions. NIST’s Handbook 44 gives 1 hectare as 10,000 square metres exactly, and 1 acre as 43,560 square feet exactly, which is 0.404 685 642 24 hectare or 4,046.856 422 4 m² (NIST). The ward of the previous program, in each unit:
| Unit | Its definition | The ward |
|---|---|---|
| Square metre | The SI unit of area | 4,001,752.60 m² |
| Hectare | 10,000 m², exactly | 400.1753 ha |
| Acre | 43,560 square feet, exactly 4,046.856 422 4 m² | 988.8546 acres |
| Bigha, katha, biswa | Set by the land records that use them | Only with the record’s definition |
Local units are another matter. A bigha, a katha or a biswa is a unit of land records, and its size is whatever the records that use it define. Haryana’s land records, for example, make a bigha 20 biswa of 20 biswansi each, a biswansi being one square karam of 57.157 inches, so that bigha is about 9,075 square feet and 4 bigha 16 biswa make an acre; for their units after settlement, the marla and the kanal, the same records use a karam of 66 inches (Haryana land records). An area in such a unit therefore means nothing without the definition it uses: convert only with the size the record itself gives, and say which definition you used.
Land Area Converter — Bigha, Guntha, Cent, Marla, Acre Convert an area between square metres, acres, hectares and the land units of a state's records; it uses a state's size only where an official source states it, and asks for your local size elsewhere.Checkpoint
Check: Land units and their definitions
Practise: measure a ward and report it in land units
Build the last two steps of a measurement report: the planar area of a ring of vertices, then the same area in hectares, acres and a local unit whose size you are given:
Code lab · 2 steps
Measure a ward and report it in land units
Measure a polygon from its vertices with the shoelace formula, then report the area in hectares, acres and a stated local unit.
Exercise · Easy · Python
Measure a polygon in a plane
Write ring_area(points) in rings.py: the area of a polygon from the (x, y) coordinates of its vertices in a plane, with the shoelace formula of the lesson's planar_area. Add x1 * y2 - x2 * y1 for every edge from one vertex to the next, the last vertex joining the first, and return half the absolute value of the sum.
Rings come both ways: some repeat the first vertex at the end and some do not, and some run clockwise, some anticlockwise. All of them must give the same area. With fewer than three distinct vertices, return 0.0.
The sample tests measure the lesson's ward, 4,000,000 m² in UTM zone 43N, a triangle and an L-shaped ring.
Starter code · rings.py
def ring_area(points):
"""The area of a polygon ring from its (x, y) vertices in a plane, whether or not the first vertex is repeated."""
# Replace this line: drop a repeated last vertex, then add x1 * y2 - x2 * y1 around the ring.
return 0.0The sample tests · test_rings.py
from rings import ring_area
WARD = [(378000, 2046000), (380000, 2046000), (380000, 2048000), (378000, 2048000)]
def test_the_ward():
"""measures the lesson's 2 km square ward in UTM metres"""
assert ring_area(WARD) == 4_000_000.0
def test_direction_and_closing():
"""gives the same area clockwise, anticlockwise and with the first vertex repeated"""
assert ring_area(WARD[::-1]) == 4_000_000.0
assert ring_area(WARD + [WARD[0]]) == 4_000_000.0
def test_other_shapes():
"""measures a right triangle and an L-shaped, concave ring"""
assert ring_area([(0, 0), (4, 0), (0, 3)]) == 6.0
assert ring_area([(0, 0), (4, 0), (4, 1), (1, 1), (1, 3), (0, 3)]) == 6.0
def test_too_few_points():
"""returns 0.0 for fewer than three distinct vertices"""
assert ring_area([]) == 0.0
assert ring_area([(1, 1), (5, 5)]) == 0.0
assert ring_area([(1, 1), (5, 5), (1, 1)]) == 0.0A hint
First drop a repeated last vertex: if len(points) > 1 and points[0] == points[-1]: points = points[:-1]. Then loop over range(len(points)) and pair each vertex with points[(i + 1) % len(points)], so the last one joins the first.
Results of the sample tests
| Test | Result | Details |
|---|
What your code printed
The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.
Exercise · Easy · Python
Report an area in hectares, acres and a local unit
Write land_units(m2, local=None) in units.py. It takes an area in square metres and returns a dict with the area in hectares, "ha", and in acres, "acre", each rounded to 4 decimals. Use the exact definitions in the file: a hectare is 10,000 m² and an acre 4,046.8564224 m².
Local units such as the bigha have no single size, so the function only converts to one when the caller says how big it is: local is a pair (name, size in m²), and the result then has a third entry under that name. Raise a ValueError when the local unit's size is not more than 0.
The sample tests use the lesson's ward, 4,001,752.6 m², the exact acre, and a local bigha defined as 1,600 square yards.
Starter code · units.py
HECTARE = 10_000.0 # square metres, exactly
ACRE = 4_046.8564224 # square metres, exactly: 43,560 international square feet
def land_units(m2, local=None):
"""An area in m² in hectares and acres, and in a local unit given as (name, size in m²), to 4 decimals."""
# Replace this line: divide by HECTARE and ACRE, and by the local unit's size when one is given.
return {"ha": 0.0, "acre": 0.0}The sample tests · test_units.py
from units import land_units
BIGHA = ("bigha", 1_600 * 0.83612736) # a local bigha defined as 1,600 square yards
def test_the_ward():
"""reports the lesson's ward, 4,001,752.6 m² on the ellipsoid, in hectares and acres"""
assert land_units(4_001_752.6) == {"ha": 400.1753, "acre": 988.8546}
def test_exact_definitions():
"""uses the exact sizes of the acre and the hectare"""
assert land_units(4_046.8564224) == {"ha": 0.4047, "acre": 1.0}
assert land_units(10_000.0)["acre"] == 2.4711
def test_a_local_unit():
"""converts to a local unit only when its size is given, as the land-area converter asks"""
assert land_units(10_000.0, BIGHA) == {"ha": 1.0, "acre": 2.4711, "bigha": 7.4749}
assert "bigha" not in land_units(10_000.0)
def test_bad_local_unit():
"""refuses a local unit without a positive size"""
try:
land_units(10_000.0, ("bigha", 0.0))
except ValueError:
return
assert False, "a size of 0 m² should raise ValueError"A hint
Build the dict with round(m2 / HECTARE, 4) and round(m2 / ACRE, 4). When local is not None, unpack it with name, size = local, check that size is positive, and add round(m2 / size, 4) under name.
Results of the sample tests
| Test | Result | Details |
|---|
What your code printed
The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.
Finally, practise the five ideas, from geodesics to land units; the session begins with the idea your answers show to be weakest:
Practice that adapts to you
Practise this lesson
12 questions about 5 ideas of this lesson. Each answer is explained at once, and what you get right is kept in this browser.
Interview questions
Warm-up (fresher): When is a planar distance in a projected CRS good enough?
When the grid’s scale factor stays close to 1 along the whole line, because a planar length is wrong by that scale minus one. In the local UTM zone the scale is within about 0.1 % of 1, so the straight line from Pune to Mumbai in zone 43N was only 2.82 m short of the 120,138.17 m geodesic. In a national conformal grid such as EPSG:7755 the same line was 1.55 % short, and in Web Mercator 5.77 % too long, which is not good enough for any measurement. So use the local UTM zone for lengths inside a city or a district, and measure geodesics on the ellipsoid for long lines, lines across zones, or anything that must match another party’s numbers.
Two questions in the module’s bank build on this: delivery distances that come out 8 % short in DuckDB, and two departments whose district areas differ by 3 %. Each lists what it tests and a first step:
Model answers open with any pass: ₹79 for 1 week ($2.99), one payment that never renews.
-
Fresher to mid
A delivery company computes trip distances in DuckDB with ST_Distance_Spheroid on points built as ST_Point(longitude, latitude). Most trips in Maharashtra come out several per cent short. Explain what happened and how you would fix and test it.
-
Mid to senior
Two departments report the area of the same district: one computed it in EPSG:7755, the other on the ellipsoid, and their figures differ by about 3 %. Which is right, how do you explain the difference, and how do you make all district statistics consistent?
Flashcards · 7 cards
Recall it, then check
Say the answer to yourself before you look: recalling it is what makes it stick.
-
What is a geodesic, and how accurate are Karney's and Vincenty's methods?
The shortest line between two points on the ellipsoid. Karney's algorithms in GeographicLib err by less than 15 nanometres; Vincenty's older formula is good to about 0.1 mm but may not converge for points almost opposite each other on the Earth.
-
How far off is a sphere, instead of the ellipsoid, for the distance from Pune to Mumbai?
The haversine formula on a sphere of radius 6,371,008.8 m gave 14.30 m more than the 120,138.17 m geodesic, 0.012 %. That is close enough for most analysis in India, but not for surveys, legal distances or checking another party's numbers.
-
How wrong is a straight-line length measured in a projected grid?
By the grid's scale along the line, minus one. Pune to Mumbai came out 2.82 m short in UTM 43N, 1.549 % short in EPSG:7755 and 5.766 % long in Web Mercator. A straight line in degrees is no distance at all.
-
Why measure areas on the ellipsoid or in an equal-area grid rather than in any projected CRS?
A conformal grid's area error is about twice its length error. Measured on the ellipsoid, a 2 km ward near Pune has 4,001,752.60 m²; its UTM zone makes it 0.044 % smaller, EPSG:7755 2.986 % smaller and Web Mercator 11.807 % larger.
-
How do QGIS, PostGIS and DuckDB decide between measuring on the ellipsoid and in the plane?
QGIS by function: $area and $length use the project's ellipsoid, area() and length() the layer's own CRS. PostGIS by type: geography on the spheroid, geometry in the units of its SRS. DuckDB by function: ST_Distance_Spheroid and ST_Area_Spheroid, which read latitude first.
-
Why did DuckDB's ST_Distance_Spheroid report 110,690 m from Pune to Mumbai instead of about 120,138 m?
The points were built as ST_Point(longitude, latitude), but the spheroid functions expect latitude first, so the axes were swapped and the answer came out 7.864 % short.
-
Which units of land area are exact, and which depend on a record's definition?
A hectare is exactly 10,000 m² and an acre exactly 43,560 square feet, or 4,046.856 422 4 m². A bigha, katha or biswa is whatever the land records using it define, so convert only with the size the record states, and say which one you used.
Key takeaways
- A distance on the Earth is a geodesic on the ellipsoid, 120,138.17 m from Pune to Mumbai; a sphere of the mean radius is 0.012 % off there, close enough for most analysis and not for surveys.
- A planar length is wrong by the grid’s scale minus one: 0.002 % in UTM 43N, −1.55 % in EPSG:7755 and +5.77 % in Web Mercator for the same line, and a length in degrees means nothing.
- Areas follow the areal scale, so measure them on the ellipsoid or in an equal-area grid: the 4,001,752.60 m² ward is 2.99 % small in EPSG:7755 and 11.81 % large in Web Mercator.
- Engines measure what the function or the type asks for:
$areaagainstarea()in QGIS,geographyagainstgeometryin PostGIS, and DuckDB’s spheroid functions, which expect latitude first. - Hectares and acres have exact definitions, 10,000 m² and 4,046.856 422 4 m²; local units such as the bigha need the definition their record uses.
Terms in this lesson
- geodesic
- The shortest line between two points on the ellipsoid; its length is the distance between them on the Earth's model, 120,138.17 m from Pune to Mumbai on WGS 84.
- haversine formula
- A formula for the great-circle distance between two points on a sphere from their latitudes and longitudes. With the Earth's mean radius it is fast and about 0.01 % off over distances such as Pune to Mumbai.
References
- C. F. F. Karney, Algorithms for geodesics, Journal of Geodesy 87 (2013) (Springer (Journal of Geodesy))
- GeographicLib: geodesics on an ellipsoid (Python documentation) (GeographicLib (C. F. F. Karney))
- GeographicLib: online geodesic calculator (GeodSolve) (GeographicLib (C. F. F. Karney))
- GeographicLib: online polygon area calculator (Planimeter) (GeographicLib (C. F. F. Karney))
- PostGIS: ST_Distance (PostGIS Project Steering Committee)
- PostGIS: ST_DistanceSphere (PostGIS Project Steering Committee)
- DuckDB spatial extension: functions (ST_Distance_Spheroid, ST_Area_Spheroid) (DuckDB Foundation)
- QGIS User Guide: expression functions ($area, area, $length, length) (QGIS Project)
- GeoPandas: GeoSeries.area (GeoPandas developers)
- NIST Handbook 44 (2026), Appendix C: general tables of units of measurement (National Institute of Standards and Technology (NIST))
- Units of measurements (WEB-HALRIS, the land records of Haryana) (Government of Haryana (jamabandi.nic.in))
Related tools
Report a problem with this lesson
Kept only in this browser. Your Learn progress