Architectural Deep-Dive: Building a High-Performance 3D Digital Twin of 500k+ Buildings in the Browser

How we engineered a lightweight, interactive 3D web map using the Google Open Buildings dataset, GEE-derived height rasters, Maplibre GL, and real-time EV charging simulations.

Introduction: The Web Digital Twin Challenge

Urban digital twins have traditionally been the domain of heavy desktop GIS packages (like ArcGIS Pro or QGIS) or massive, proprietary enterprise platforms. These setups require high-end workstations and proprietary licenses, creating a barrier to entry for urban planners, municipal bodies, and citizens in developing urban centers.

Building a fully interactive, web-based digital twin that runs at 60 FPS in a standard mobile or desktop browser is a daunting engineering challenge. To visualize a city the size of Chennai (covering over 1,189 square kilometers of the Chennai Metropolitan Area), we had to process and render over 525,000 buildings, overlay spatial datasets, and support real-time simulations.

A naive representation of this data in standard GeoJSON format exceeds 1.5 GB—a size that would crash any web browser, blow past mobile data limits, and stall CPU main threads.

In this article, we details the end-to-end engineering pipeline we built to solve this problem:

Ingestion & Streaming Filter: Streaming and filtering millions of footprints on-the-fly.
Vectorized Spatial Joins: Merging building shapes with satellite height data using matrix math.
Urban Semantics: Running a spatial rule-engine to classify building use.
Data Compacting: Compressing the payload by 96% (from 1.5 GB down to 60 MB).
High-Performance Web 3D Rendering: Delivering a smooth, hardware-accelerated client app featuring real-time simulators and custom night maps.

The System Architecture

The following diagram illustrates the architecture of our geospatial processing and visualization pipeline:

Figure 1 • 3D Digital Twin Architecture (525k Buildings)

Interactive 5-Stage Spatial Data Engineering Pipeline

1. Streaming & Spatial Filtering

Bypasses heavy downloads by filtering candidate footprints on-the-fly.

Python / Shell Pipe

Raw Dataset Google Open Buildings CSV.gz

Stream Command curl | gunzip | awk bounding filter

Index Engine Shapely Prepared Geometries (prep)

CMA Candidate Count 525,000+ Footprints

Engineering Context: Command-line pipe filters coordinate streams before decompressing into Python memory. Prepared geometries accelerate complex polygon containment checks.

1. Streaming & Spatial Filtering of Big Geospatial Data

The Google Open Buildings dataset is an open-source resource containing millions of AI-detected building footprints. However, download sizes are immense. The tile covering Chennai and surrounding parts of Southern India (3a5_buildings.csv.gz) is a gzipped CSV containing millions of building rows.

To avoid downloading gigabytes of unnecessary data onto local machines, we built a streaming pipeline in Python that filters and processes coordinates on the fly.

On-the-Fly Stream Filtering

We combined shell utilities (curl, gunzip, and awk) directly within a subprocess pipe. This allowed us to filter rows by their bounding boxes before they were even decompressed into python memory:

stream_filtering.py PYTHON

cmd = (

f"curl -L -s '{FOOTPRINT_URL}' | "

f"gunzip -c | "

f"awk -F, 'NR==1 || ($1 >= {MIN_LAT} && $1 <= {MAX_LAT} && $2 >= {MIN_LON} && $2 <= {MAX_LON})' "

f"> {FILTERED_CSV_PATH}"

)

curl pulls the compressed stream.
gunzip decompresses the stream line-by-line.
awk filters coordinates based on bounding box limits.
Result: We shrunk a massive tile down to a temporary local CSV containing only the candidate structures inside the rectangular boundary of Chennai.

Boundary Intersection via Prepared Geometries

A bounding box is a rectangular approximation. The actual political boundary of the Chennai Metropolitan Area (CMA) is a highly complex, irregular polygon. Running standard point-in-polygon checks (like polygon.contains(point)) for hundreds of thousands of candidate buildings is computationally expensive.

To optimize this, we utilized Prepared Geometries from the shapely library:

prepared_geometry_containment.py PYTHON

from shapely.geometry import shape, Point

from shapely.prepared import prep

# Load political boundary

with open("CMA.geojson", 'r') as f:

cma_geom = shape(json.load(f)["features"][0]["geometry"])

# Compile boundary into an optimized prepared geometry

prepared_cma = prep(cma_geom)

# Centroid check inside loop

pt = Point(lon, lat)

if prepared_cma.contains(pt):

# Keep footprint

shapely.prepared.prep sets up an in-memory spatial index (using a modified R-tree structure) that accelerates geometric containment queries. By doing a quick containment check on the building's centroid before parsing its complex polygon structure (via WKT parser), we completed the filtering of 525,000+ exact footprints in seconds.

2. Vectorized Spatial Joins: Merging Footprints and Raster Heights

The Google Open Buildings dataset provides footprints, but to create a true 3D Digital Twin, we needed heights. We used a building-height raster generated from Google Earth Engine (GEE).

Raster data is stored as a 2D grid of pixels (each pixel corresponding to a real-world cell, in our case, representing height in meters). Vector data (the building footprints) consists of latitude/longitude coordinates. Intersecting half a million vector shapes with a high-resolution GeoTIFF raster using traditional spatial libraries is slow.

Slicing with Vectorized Coordinates

Instead of looping and sampling the raster pixel-by-pixel, we vectorized the coordinate transform using numpy and rasterio:

vectorized_height_intersection.py PYTHON

import numpy as np

import rasterio

from rasterio.transform import rowcol

# Read raster band into a memory array

with rasterio.open("cma_building_height_2023.tif") as src:

band1 = src.read(1)

height_dim, width_dim = band1.shape

# Extract lons and lats from all features

lons = [feat["properties"]["centroid_lon"] for feat in features]

lats = [feat["properties"]["centroid_lat"] for feat in features]

# Calculate pixel indices for all centroids simultaneously using matrix math

rows, cols = rowcol(src.transform, lons, lats)

rows = np.array(rows)

cols = np.array(cols)

# Boundary mask to ensure coordinates fall inside the raster bounds

valid_mask = (rows >= 0) & (rows < height_dim) & (cols >= 0) & (cols < width_dim)

# Sample height values using NumPy indexing

sampled_heights = np.full(len(features), 3.0, dtype=np.float32)

sampled_heights[valid_mask] = band1[rows[valid_mask], cols[valid_mask]]

Rather than querying the filesystem or raster structure iteratively, rowcol uses matrix algebra on the affine transformation to map coordinates to pixel rows/columns in a single execution. The sampling step is a direct array index band1[rows, cols], leveraging optimized C-underpinnings of NumPy.

Any invalid pixels (e.g. NaNs, nodata values, or heights below 0) were normalized to a default single-story height of 3.0m.

Figure 2 • Urban Heuristic Rule-Tree Simulator

Test real-time building floor count & UDPFI classification

Building Height 18.5m (5 Floors)

3.0m (1 Flr) 21m (6 Flrs) 45m (13 Flrs)

Footprint Area 450 m²

50 m² 600 m² 1500 m²

UDPFI Urban Land-Use Class

[AM] Apartments (Medium Rise)

Compact GeoJSON Output (96% Compressed)

{
  "type": "Feature",
  "properties": {
    "h": 18.5,
    "a": 450,
    "u": "AM"
  }
}

Short keys (h, a, u) and 2-character land-use codes allow Maplibre GL vertex shaders to style 525,000 extruded 3D features dynamically in WebGL memory.

3. Classifying Building Use via Spatial Heuristics

Real-world digital twins must map building functions (Residential, Commercial, Industrial) to help urban planners. Since open building datasets rarely include land-use labels, we engineered a spatial heuristic classifier in Python.

Our rule-engine utilizes two spatial attributes: building height (representing capacity) and footprint area (representing floorplate size). First, we estimate the floor count:

$\text{Floors} = \max\left(1, \text{round}\left(\frac{\text{Height}}{3.5\text{m}}\right)\right)$

Using these values, we ran a heuristic tree aligned with local urban typologies:

heuristic_use_classifier.py PYTHON

def classify_building(height, area):

floors = max(1, round(height / 3.5))

if floors == 1:

if area < 150:

return "Residential (Low Density)"

elif area < 600:

return "Commercial / Retail"

else:

return "Industrial / Warehouse"

elif floors <= 3:

if area < 250:

return "Residential / Independent House"

elif area < 600:

return "Apartments / Mixed-Use"

else:

return "Commercial / Office / Retail"

elif floors <= 6:

if area < 500:

return "Apartments (Medium Rise)"

elif area < 1200:

return "Commercial / Office"

else:

return "Institutional / Public Building"

else: # 7+ floors (High Rise)

if area < 800:

return "Apartments (High Rise)"

else:

return "Commercial Office Tower / Corporate Hub"

This classification translates building geometries into urban semantics, allowing the frontend to style the city according to UDPFI (Urban Development Plans Formulation and Implementation) town planning color standards.

4. Shrinking 1.5 GB GeoJSON to 60 MB

Loading a 1.5 GB GeoJSON in a browser is impractical. We implemented a multi-stage optimization pipeline to compress the geospatial payload by over 96%:

1. Coordinate Precision Reduction

Double-precision coordinates in GeoJSON (e.g., 80.146249102948123) contain up to 15 decimal places. That represents sub-millimeter precision—far beyond the resolution of satellite datasets.

6 decimal places (80.146249) provides 10 cm accuracy, which is perfect for building outlines.
By rounding all coordinates to 6 decimal places, we cut millions of characters:

round_coords.py PYTHON

def round_coords(coords):

if isinstance(coords[0], list):

return [round_coords(c) for c in coords]

return [round(coords[0], 6), round(coords[1], 6)]

2. Attribute Compacting (Short-Coding)

GeoJSON properties add redundant string overhead. Key names like height and estimated_use repeat in every single building feature. We compacted keys and mapped values to two-character codes:

height $\rightarrow$ h (rounded to 1 decimal place)
area $\rightarrow$ a (converted to integer)
estimated_use $\rightarrow$ u (using codes like RL for Residential Low, CR for Commercial Retail, IW for Industrial)

attribute_compacting.json JSON

/* BEFORE: 215 bytes */

{

"type": "Feature",

"properties": {

"height": 18.42319,

"area": 420.1582,

"estimated_use": "Apartments / Mixed-Use"

"geometry": { ... }

}

/* AFTER: 98 bytes */

{

"type": "Feature",

"properties": {

"h": 18.4,

"a": 420,

"u": "AM"

"geometry": { ... }

}

3. Minification and Gzip Level 9

We dumped the JSON without any indentation or whitespaces using python's separators (’,’, ’:’), and then compressed it using gzip's maximum compression level (compresslevel=9):

geojson_gzip_export.py PYTHON

with open(temp_output_path, 'w', encoding='utf-8') as f:

json.dump(compact_data, f, separators=(',', ':'))

Result: The file size dropped from a raw 1.47 GB to a 60 MB .geojson.gz archive.

Figure 3 • Non-Linear EV Battery Charging Simulator

Calculates BMS cell protection taper phases (>80% State of Charge)

Vehicle Profile

Charger Type

Initial State of Charge (SoC) 20%

Target State of Charge (SoC) 90%

Estimated Charging Time

59m total duration

Fast Phase (20% → 80%): 32 mins

BMS Taper Phase (80% → 90%): 27 mins (5x Slowdown)

BMS Physics Note: When charging above 80% on DC Fast Chargers, the Battery Management System (BMS) tapers current to prevent lithium plating and thermal degradation.

5. High-Performance Client-Side Rendering & EV Simulation

In the client browser, we integrated a combination of technologies to load, decompress, and render the digital twin at high framerates.

On-the-Fly Client-Side Decompression

Rather than downloading an uncompressed file or relying on server-side middleware (which isn't always supported on basic CDNs), the web app fetches the .geojson.gz array buffer and decompresses it client-side using fflate:

client_fflate_decompression.js JAVASCRIPT

fetch('cma_buildings_3d.geojson.gz')

.then(res => res.arrayBuffer())

.then(buf => {

// High-speed sync decompression in JS

const decompressed = fflate.gunzipSync(new Uint8Array(buf));

const jsonText = new TextDecoder().decode(decompressed);

const data = JSON.parse(jsonText);

// Feed directly to Maplibre

map.addSource('buildings', { type: 'geojson', data: data });

});

fflate is a lightweight, pure-JS decompression library that is significantly faster than standard browser-based decoding libraries. It completes decompression and JSON parsing of half a million features in less than 2 seconds on modern laptops.

WebGL-Accelerated 3D Rendering

For rendering, we used Maplibre GL JS, an open-source WebGL/WebGPU-based library. We defined a fill-extrusion layer, mapping building heights (h) and use classifications (u) to WebGL paint properties:

maplibre_3d_extrusion.js JAVASCRIPT

map.addLayer({

'id': 'buildings-3d',

'type': 'fill-extrusion',

'source': 'buildings',

'paint': {

// Dynamic extrusion based on property height

'fill-extrusion-height': ['get', 'h'],

'fill-extrusion-base': 0,

// Color mapping using Maplibre expression evaluation

'fill-extrusion-color': [

'match', ['get', 'u'],

'RL', '#facc15', // Residential Low (Yellow)

'CO', '#2563eb', // Commercial Office (Blue)

'IW', '#a855f7', // Industrial (Violet)

'IP', '#ef4444', // Institutional (Red)

'#94a3b8' // Default

'fill-extrusion-opacity': 0.85

}

});

By delegating geometry extrusion and styling to vertex shaders on the GPU, the browser handles hundreds of thousands of polygons without burdening the JavaScript execution stack.

Dynamic Night Mode Paint Engine

To prevent map flashes when toggling between themes, we implemented a paint engine that changes the colors of the underlying Maplibre vector tiles on the fly, instead of changing the map style:

night_mode_paint_engine.js JAVASCRIPT

function applyBaseThemePaint() {

const isNight = activeTheme === 'night';

// Change background style

map.setPaintProperty('background', 'background-color', isNight ? '#000000' : '#0b0f19');

// Convert water features into glowing cyan

map.setPaintProperty('water', 'fill-color', isNight ? '#00e5ff' : '#0284c7');

map.setPaintProperty('water', 'fill-opacity', isNight ? 0.35 : 0.6);

// Change road casings to fit dark neon style

map.setPaintProperty('road_primary', 'line-color', isNight ? '#2a2a2a' : '#3c4858');

}

Non-Linear EV Battery Charging Simulator

We overlaid the Digital Twin with India-wide EV charging locations. To show how a digital twin can run interactive models, we built a client-side EV charging duration simulator.

Lithium-ion batteries do not charge linearly. Fast charging is effective up to 80% State of Charge (SoC). Beyond 80%, the battery management system (BMS) reduces charging speeds to prevent cell degradation.

Our simulator calculates charging duration by splitting the charging profile into fast and slow phases when the target charge exceeds 80%:

ev_battery_charging_simulator.js JAVASCRIPT

const curPct = parseInt(currentBatterySlider.value); // e.g., 20%

const tgtPct = parseInt(targetBatterySlider.value); // e.g., 90%

const capacity = activeBatteryCapacity; // e.g., 45 kWh (Tata Curvv.ev)

const power = chargerPower; // e.g., 50 kW (DC Fast Charger)

let chargeHours = 0;

if (isDC && tgtPct > 80) {

// 1. Fast Charging Phase (up to 80%)

const fastCapacity = capacity * (Math.min(80, tgtPct) - curPct) / 100;

chargeHours += fastCapacity > 0 ? (fastCapacity / power) : 0;

// 2. Slow Charging Phase (80% to target)

const slowCapacity = capacity * (tgtPct - Math.max(80, curPct)) / 100;

// Charge speed drops to 20% of maximum power (5x slowdown)

chargeHours += slowCapacity > 0 ? (slowCapacity / (power * 0.2)) : 0;

} else {

// Linear calculation for AC charging or sub-80% DC charging

const capacityNeeded = capacity * (tgtPct - curPct) / 100;

chargeHours = capacityNeeded / power;

}

When users hover over a charging station, the map calculates exact charging times based on the selected vehicle's battery capacity, the station's charger type (AC Type-2 vs DC CCS2), and the user's targeted charge levels.

Conclusion: Key Takeaways for Geospatial Engineers

Our project demonstrates that complex 3D digital twins can run effectively on standard web browsers. When building web-based geospatial tools, consider these key strategies:

Optimize at Source: Stream and filter datasets using command-line pipelines (like awk or grep) before loading them into memory.
Vectorize Spatial Operations: Avoid loops in python. Use libraries like NumPy and Rasterio to perform affine transformations on arrays.
Minimize Payload Size: Round coordinate decimals, shorten keys, map variables to small codes, and use high-level Gzip compression.
Use GPU Acceleration: Delegate 3D rendering to WebGL/WebGPU by utilizing Maplibre GL.
Add Interactive Simulations: Use client-side logic to run models, keeping your web applications fast and engaging.

Published in Research & Policy Essays

Discuss this methodology →