Address to site: parcel area, distances and shares

The longitude and latitude box that holds a radius in meters around a point.

Δlat = r ÷ 111320; Δlon = r ÷ (111320 × max(0.01, cos lat))

64-bit, with the typed figure 111320 m per degree for both directions and the cosine floored at 0.01. It is a search window, so the approximation only widens or narrows what is looked at.

export function lonLatBox(lat, lon, radiusM) {
  const dLat = radiusM / 111320, dLon = radiusM / (111320 * Math.max(0.01, Math.cos((lat * Math.PI) / 180)));
  return { minX: lon - dLon, minY: lat - dLat, maxX: lon + dLon, maxY: lat + dLat };
}

The signed area of one polygon ring.

A = ½ Σ (xᵢ·yᵢ₊₁ − xᵢ₊₁·yᵢ)

The shoelace sum on the coordinates as read, with no shift to a local origin. Each term subtracts two products of the size of the coordinates squared, while the area is far smaller, so leading digits cancel. For a parcel about 30 m across given in longitude and latitude near 122° and 38°, the products are about 3 × 10¹⁰ times twice the area, which leaves about 5 significant digits.

const ringArea = (r) => { let a = 0; for (let i = 0, n = r.length; i < n; i++) { const p = r[i], q = r[(i + 1) % n]; a += p[0] * q[1] - q[0] * p[1]; } return a / 2; };

The area of a polygon with holes, scaled to square meters.

A = (Σ |outer ring| − Σ |holes|) × kx × ky

Ring areas are added and subtracted in 64-bit and multiplied by the two meters-per-unit factors. In longitude and latitude those are 111320 × cos(latitude of the query point) and 111320, a flat approximation applied to the whole polygon.

function geomArea(g, kx, ky) { let a = 0; for (const poly of g.c) { a += Math.abs(ringArea(poly[0])); for (const h of poly.slice(1)) a -= Math.abs(ringArea(h)); } return a * kx * ky; }

The factors that turn longitude and latitude differences into meters, and the search window.

kx = 111320 × cos(lat); ky = 111320

One cosine at the query point's latitude in 64-bit. 111320 is a typed spherical figure; the length of a degree on the ellipsoid differs from it by up to about 0.7 % depending on latitude and direction, and the code does not correct for that.

const kx = geographic ? 111320 * Math.cos((y * Math.PI) / 180) : 1, ky = geographic ? 111320 : 1;
const r = geographic ? { dx: radius / kx, dy: radius / ky } : { dx: radius, dy: radius };
const area = { minX: x - r.dx, minY: y - r.dy, maxX: x + r.dx, maxY: y + r.dy };

The distance from the point to each building and the building's footprint area.

d = √(((cx − x)·kx)² + ((cy − y)·ky)²), c = center of the building's bounding box

Math.hypot of the scaled differences. The distance is to the center of the building's bounding box, not to its nearest wall or its true centroid.

for (const { f, over } of feats('buildings')) {
  const c = centroid(f.g), d = Math.hypot((c.x - x) * kx, (c.y - y) * ky);
  if (d > radius) continue;
  const rec = mapFeature('buildings', f.p, over), levels = num(rec.levels);
  site.buildings.push({ record: { height: num(rec.height), levels, kind: text(rec.kind), name: text(rec.name) }, area: geomArea(f.g, kx, ky), d, x: c.x, y: c.y });
}
site.buildings.sort((a, b) => a.d - b.d);

The share of the parcel that lies in each zoning district, soil unit and flood zone.

share = points in the zone ÷ points in the parcel, on a 40 × 40 grid over the parcel's bounding box

A count of sample points, not an area of overlap, so a share can only take steps of one over the number of points that fall inside the parcel. Each point takes the first polygon that contains it. With no parcel, the single query point is the only sample.

const samples = [];
if (site.parcel) {
  const b = site.parcel.bounds, nx = grid, ny = grid;
  for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) { const px = b.minX + ((i + 0.5) * (b.maxX - b.minX)) / nx, py = b.minY + ((j + 0.5) * (b.maxY - b.minY)) / ny; if (best.f.g.c.some((poly) => inPoly(px, py, poly))) samples.push([px, py]); }
}
if (!samples.length) samples.push([x, y]);
const share = (kind, key) => {
  const list = feats(kind), tally = new Map();
  for (const [px, py] of samples) { const hit = list.find(({ f }) => inGeom(px, py, f.g)); const k = hit ? key(hit) : null; if (k !== null) tally.set(k.id, { ...k, n: (tally.get(k.id)?.n || 0) + 1 }); }
  return [...tally.values()].sort((a, b) => b.n - a.n).map((v) => ({ ...v, share: v.n / samples.length }));
};