DEV Community

jj yang
jj yang

Posted on AI-assisted

Finding Nearby Cities Along a Polyline in TypeScript

Finding Nearby Cities Along a Polyline in TypeScript

Most “nearby” features start with a point.

This one starts with a path.

Given an anchor point and a polyline made of latitude/longitude coordinates, we want to return a small list of real cities that satisfy both conditions:

  1. The city is close to the anchor point.
  2. The city is close to the actual line, not merely to one sampled vertex.

I ran into this while building map tooling for AstroCarto. The product context is an interactive astrocartography map, but the engineering problem is generic. The same pattern works for delivery corridors, hiking routes, railway lines, service areas, and coverage maps.

Define “nearby” as two constraints

A city can be close to the location a user is exploring but far from the selected line. It can also be close to the line while being too far from the current area of interest.

So the query is:

distance(city, anchor) <= anchorRadius
AND
distance(city, polyline) <= lineRadius
Enter fullscreen mode Exit fullscreen mode

The result should be allowed to contain zero, one, two, or several cities. Returning a distant city just to fill a UI slot creates false precision.

A useful empty result is better than a complete-looking but misleading list.

Why not scan everything or query a database?

A full in-memory scan is easy to write, and PostGIS is an excellent choice for dynamic or very large datasets. Neither is automatically wrong.

The static-client approach is useful when:

  • the source data changes infrequently;
  • the query is local and bounded;
  • the browser can load only the relevant region;
  • the final calculation is small enough for the device;
  • a map interaction should not depend on a server round trip.

For changing data, millions of records, polygons, joins, or access-controlled results, a spatial database is usually the better boundary. This article is about the smaller read-heavy case.

A minimal, runnable TypeScript implementation

The following example deliberately uses synthetic data and illustrative thresholds. They are not live configuration values.

It uses a spherical point-to-segment calculation rather than a flat x/y approximation.

Save it as nearby-polyline-example.ts and run it with:

npx tsx nearby-polyline-example.ts
Enter fullscreen mode Exit fullscreen mode
type Point = { lat: number; lon: number };
type City = Point & { id: number; name: string; population: number };
type Vector = readonly [number, number, number];

const EARTH_RADIUS_KM = 6_371.0088;
const DEG = Math.PI / 180;

const clamp = (value: number, min: number, max: number) =>
  Math.max(min, Math.min(max, value));

function toVector(point: Point): Vector {
  const lat = point.lat * DEG;
  const lon = point.lon * DEG;
  const cosLat = Math.cos(lat);

  return [
    cosLat * Math.cos(lon),
    cosLat * Math.sin(lon),
    Math.sin(lat),
  ];
}

function dot(a: Vector, b: Vector) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

function cross(a: Vector, b: Vector): Vector {
  return [
    a[1] * b[2] - a[2] * b[1],
    a[2] * b[0] - a[0] * b[2],
    a[0] * b[1] - a[1] * b[0],
  ];
}

function normalize(value: Vector): Vector | null {
  const length = Math.hypot(...value);

  return length < 1e-15
    ? null
    : [value[0] / length, value[1] / length, value[2] / length];
}

function angle(a: Vector, b: Vector) {
  return Math.acos(clamp(dot(a, b), -1, 1));
}

function haversineKm(a: Point, b: Point) {
  const dLat = (b.lat - a.lat) * DEG;
  const dLon = (b.lon - a.lon) * DEG;
  const lat1 = a.lat * DEG;
  const lat2 = b.lat * DEG;

  const h =
    Math.sin(dLat / 2) ** 2 +
    Math.cos(lat1) *
      Math.cos(lat2) *
      Math.sin(dLon / 2) ** 2;

  return 2 * EARTH_RADIUS_KM * Math.asin(Math.sqrt(clamp(h, 0, 1)));
}

function pointToSegmentKm(point: Point, start: Point, end: Point) {
  const p = toVector(point);
  const a = toVector(start);
  const b = toVector(end);
  const segmentAngle = angle(a, b);

  if (segmentAngle < 1e-15) {
    return angle(p, a) * EARTH_RADIUS_KM;
  }

  const normal = normalize(cross(a, b));

  if (!normal) {
    return Math.min(angle(p, a), angle(p, b)) * EARTH_RADIUS_KM;
  }

  // Project p onto the great-circle plane.
  const projection = normalize([
    p[0] - dot(p, normal) * normal[0],
    p[1] - dot(p, normal) * normal[1],
    p[2] - dot(p, normal) * normal[2],
  ]);

  let best = Math.min(angle(p, a), angle(p, b));

  for (const candidate of projection
    ? [projection, [-projection[0], -projection[1], -projection[2]] as const]
    : []) {
    // Keep the projection only if it lies on the minor arc.
    const onMinorArc =
      Math.abs(angle(a, candidate) + angle(candidate, b) - segmentAngle) <
      1e-8;

    if (onMinorArc) {
      best = Math.min(best, angle(p, candidate));
    }
  }

  return best * EARTH_RADIUS_KM;
}

function pointToPolylineKm(point: Point, line: Point[]) {
  let best = Number.POSITIVE_INFINITY;

  for (let i = 0; i < line.length - 1; i += 1) {
    // Split rendering seams at the antimeridian before calling this function.
    best = Math.min(best, pointToSegmentKm(point, line[i], line[i + 1]));
  }

  return best;
}

function findNearbyCities(
  anchor: Point,
  line: Point[],
  cities: City[],
  // Illustrative example values only.
  options = { maxAnchorKm: 300, maxLineKm: 100 },
  limit = 5,
) {
  return cities
    .map((city) => ({
      ...city,
      anchorKm: haversineKm(anchor, city),
      lineKm: pointToPolylineKm(city, line),
    }))
    .filter(
      (city) =>
        city.anchorKm <= options.maxAnchorKm &&
        city.lineKm <= options.maxLineKm,
    )
    .sort(
      (a, b) =>
        a.lineKm - b.lineKm ||
        b.population - a.population ||
        a.id - b.id,
    )
    .slice(0, limit)
    .map(({ id, name, anchorKm, lineKm }) => ({
      id,
      name,
      anchorKm: Math.round(anchorKm),
      lineKm: Math.round(lineKm),
    }));
}

const anchor = { lat: 51.5074, lon: -0.1278 };

const line = [
  { lat: 50.8, lon: -1.8 },
  { lat: 51.5074, lon: -0.1278 },
  { lat: 52.6, lon: -1.7 },
];

const cities: City[] = [
  {
    id: 1,
    name: "London",
    lat: 51.5074,
    lon: -0.1278,
    population: 9_000_000,
  },
  {
    id: 2,
    name: "Southampton",
    lat: 50.9097,
    lon: -1.4044,
    population: 250_000,
  },
  {
    id: 3,
    name: "Birmingham",
    lat: 52.4862,
    lon: -1.8904,
    population: 1_100_000,
  },
  {
    id: 4,
    name: "Manchester",
    lat: 53.4808,
    lon: -2.2426,
    population: 550_000,
  },
];

console.log(findNearbyCities(anchor, line, cities));
Enter fullscreen mode Exit fullscreen mode

Output:

[
  { id: 1, name: 'London', anchorKm: 0, lineKm: 0 },
  { id: 2, name: 'Southampton', anchorKm: 111, lineKm: 6 },
  { id: 3, name: 'Birmingham', anchorKm: 162, lineKm: 18 }
]
Enter fullscreen mode Exit fullscreen mode

The code is intentionally small. It does not include a tile index, a cache, or a metropolitan clustering policy. Those are delivery and product decisions layered on top of the geometry.

The spherical geometry

Latitude and longitude are angular coordinates on a sphere. A flat distance formula can be acceptable for a very small local map, but it becomes misleading at high latitudes and around the date line.

For a coordinate (latitude, longitude), convert it to a unit vector:

v = (
  cos(latitude) * cos(longitude),
  cos(latitude) * sin(longitude),
  sin(latitude)
)
Enter fullscreen mode Exit fullscreen mode

For a segment with vectors a and b, the great-circle plane has a normal:

n = normalize(a × b)
Enter fullscreen mode Exit fullscreen mode

Project the candidate point p onto that plane:

q = normalize(p - (p · n)n)
Enter fullscreen mode Exit fullscreen mode

The projection is valid only when it falls on the finite minor arc. The implementation checks:

angle(a, q) + angle(q, b) ≈ angle(a, b)
Enter fullscreen mode Exit fullscreen mode

If the projection is outside the segment, the closest point is one of the endpoints.

For a polyline, repeat this for every segment and keep the minimum distance.

If the original path is an analytic curve rather than a geodesic, sample it densely enough. The accuracy of the result cannot exceed the accuracy of the sampled line.

Deliver a large city corpus without shipping it all at once

The geometry above assumes that all cities are already in memory. A global map usually needs a better delivery strategy.

The general pipeline is:

versioned manifest
  -> choose assets near the current area
  -> fetch only those assets
  -> reuse HTTP and memory caches
  -> filter by anchor distance
  -> calculate exact spherical line distance
  -> rank and cluster
  -> return a small result set
Enter fullscreen mode Exit fullscreen mode

The data build can happen offline:

  1. Download a public city source such as GeoNames.
  2. Validate coordinates, IDs, and country codes.
  3. Remove records that are not useful as independent destinations.
  4. Assign a simple, explainable importance field.
  5. Partition the result into geographic assets.
  6. Publish a versioned manifest and compressed files.

The browser only needs a compact record such as:

type CityRecord = {
  id: number;
  name: string;
  countryCode: string;
  latitude: number;
  longitude: number;
  importance: number;
};
Enter fullscreen mode Exit fullscreen mode

Keep the source version explicit. A mutable latest asset makes debugging and rollback much harder.

For a local interaction, select assets near the anchor first. If the first pass produces too few results, perform one wider pass and reuse the data already loaded. For a long route, select assets around the line envelope rather than around one anchor point.

The important optimization is not a clever sort. It is reducing the number of candidate records before doing exact geometry.

Ranking: closest is not always most useful

Distance alone can fill the result with several suburbs from the same metropolitan area.

A more useful ranking separates:

  1. Distance to the current anchor.
  2. Distance to the selected line.
  3. City representativeness.
  4. Regional diversity.

One deterministic ordering might be:

proximity bucket
-> importance
-> anchor distance
-> line distance
-> stable city ID
Enter fullscreen mode Exit fullscreen mode

After sorting, merge cities that fall within an application-defined metropolitan radius and keep one representative. The radius should be validated against the geography of the target market; it should not be presented as a universal truth.

Always allow the result to contain fewer cities than the UI limit. Sparse geography is information.

Edge cases worth testing

Coordinate order

Choose one internal convention and convert at the boundary. GeoJSON uses [longitude, latitude], while many application models use [latitude, longitude]. A swapped pair can produce plausible-looking but completely wrong results.

Antimeridian crossings

A segment from 179°E to 179°W is a short segment near the date line, not a trip around the entire world. Split rendering seams at ±180° or explicitly unwrap longitudes before calculating distances.

Degenerate segments

Identical endpoints do not define a useful great-circle normal. Fall back to point distance. Nearly antipodal endpoints also need a numerical fallback.

Stale asynchronous results

When a user selects another line while assets are loading, associate each request with a request ID or an AbortController. Apply a result only if it still belongs to the current selection.

Asset failures

An optional city asset should not make the map unusable. Keep the selected line and its basic information visible, show a retry action, and remove failed promises from the cache so retrying is possible.

Empty regions

Oceans, polar regions, and sparsely populated areas can legitimately return no cities. Do not silently substitute a distant city just to make the interface look complete.

How to benchmark it honestly

Measure the stages separately:

  1. Network transfer.
  2. JSON parsing.
  3. Candidate filtering.
  4. Exact line-distance calculation.
  5. Ranking and rendering.

In my local evaluation, fixtures covered dense cities, sparse regions, polar areas, and date-line geometry. The in-memory filtering and geometry stage stayed within a small interactive budget after assets had loaded.

That result is useful for comparing algorithm changes, but it is not a mobile performance guarantee. Browser network conditions, device CPU, map rendering, and line sampling density should be measured separately.

When a spatial database is the better choice

Use PostGIS or another server-side spatial index when:

  • records change frequently;
  • the dataset is much larger than a browser should hold;
  • queries involve polygons, buffers, joins, or permissions;
  • multiple clients need a consistent view of the data;
  • you need server-side observability and access control.

For a larger static archive, vector tiles, PMTiles, or FlatGeobuf may be a better delivery format. H3 and geohash indexes are useful for coarse candidate selection, but the final distance to the actual line still needs a proper geometry calculation.

The product lesson

The original requirement came from a map interaction where a selected line needed a few geographically meaningful places to inspect next. The domain was specific, but the engineering lesson was not:

bounded data
+ local asset loading
+ spherical geometry
+ deterministic ranking
+ graceful empty and error states
Enter fullscreen mode Exit fullscreen mode

That combination is often enough for a responsive browser-side map feature without turning every interaction into a database problem.

If you want to see the broader interactive map product that motivated this work, you can explore AstroCarto.

Disclosure: AstroCarto was the product context for this build note. The algorithm is intended to be reusable outside astrology.

City data examples use public geographic information from GeoNames; retain the applicable CC BY 4.0 attribution when using that source.

Top comments (0)