Speed up address and road diffs from hours to seconds, gzip GeoJSON downloads, and render map points on canvas with viewport-only OSM/county layers

compare-addresses.py

@@ -99,9 +99,12 @@ class AddressComparator:
         Returns the number of pairs found."""
         if not new_addresses or not removed_addresses:
             return 0
+        # Index the unmatched OSM side by house number; a pair needs the same one.
         geoms = [a['geometry'] for a in removed_addresses]
-        tree = STRtree(geoms)
+        by_num: Dict[str, List[int]] = {}
-        radius_deg = radius_meters / 111000.0
+        for j, osm in enumerate(removed_addresses):
+            by_num.setdefault(self._normalize_multi(osm.get('addr:housenumber')), []).append(j)
+        street_cache: Dict[int, str] = {}
         paired = set()
         count = 0
         for addr in new_addresses:
@@ -111,13 +114,13 @@ class AddressComparator:
             if not num or not street:
                 continue
             best, best_dist = None, float('inf')
-            for j in tree.query(addr['geometry'].buffer(radius_deg)):
+            for j in by_num.get(num, []):
                 if j in paired:
                     continue
                 osm = removed_addresses[j]
-                if self._normalize_multi(osm.get('addr:housenumber')) != num:
+                if j not in street_cache:
-                    continue
+                    street_cache[j] = self._normalize_street_name(str(osm.get('addr:street', '')))
-                osm_street = self._normalize_street_name(str(osm.get('addr:street', '')))
+                osm_street = street_cache[j]
                 if not osm_street or osm_street == street:
                     continue
                 osm_unit = self._normalize_multi(osm.get('addr:unit'))
@@ -162,97 +165,83 @@ class AddressComparator:
         if local_gdf.crs != osm_gdf.crs:
             osm_gdf = osm_gdf.to_crs(local_gdf.crs)
 
-        osm_index = STRtree(osm_gdf.geometry.tolist())
+        # Normalise every address once and index OSM by (number, street):
+        # a county address can only ever match an OSM address with the same
+        # key, so there's no need to scan everything within the tolerance
+        # (thousands of points in dense areas like The Villages).
+        def keys(gdf):
+            nums = [self._normalize_multi(v) for v in gdf.get('addr:housenumber', pd.Series([''] * len(gdf)))]
+            streets = [self._normalize_street_name(str(v)) for v in gdf.get('addr:street', pd.Series([''] * len(gdf)))]
+            units = [self._normalize_multi(v) for v in gdf.get('addr:unit', pd.Series([None] * len(gdf)))]
+            return nums, streets, units
+
+        osm_nums, osm_streets, osm_units = keys(osm_gdf)
+        local_nums, local_streets, local_units = keys(local_gdf)
+        osm_geoms = osm_gdf.geometry.tolist()
+        local_geoms = local_gdf.geometry.tolist()
+        osm_props = osm_gdf.drop(columns='geometry').to_dict('records')
+        local_props = local_gdf.drop(columns='geometry').to_dict('records')
+
+        osm_by_key: Dict[Tuple[str, str], List[int]] = {}
+        for i, (num, street) in enumerate(zip(osm_nums, osm_streets)):
+            if street:
+                osm_by_key.setdefault((num, street), []).append(i)
 
         existing_addresses = []
         new_addresses = []
+        matched_osm = set()
 
-        for _, local_row in local_gdf.iterrows():
+        for i, local_point in enumerate(local_geoms):
-            local_point = local_row.geometry
+            street = local_streets[i]
-            nearby = osm_index.query(local_point.buffer(self.tolerance_deg))
+            candidates = osm_by_key.get((local_nums[i], street), []) if street else []
+            local_unit = local_units[i]
 
-            best_match = None
+            # Existing: nearest OSM address with the same unit (or both
-            min_distance = float('inf')
+            # without one) within the tolerance.
+            best_match, min_distance = None, float('inf')
+            # Matched (not potentially removed): the same, except an OSM
+            # address with no unit also covers a county address with one.
+            nearest_cover, cover_distance = None, float('inf')
+            for j in candidates:
+                distance = local_point.distance(osm_geoms[j])
+                if distance * 111000.0 > self.tolerance_meters:
+                    continue
+                osm_unit = osm_units[j]
+                if (local_unit == osm_unit) and distance < min_distance:
+                    best_match, min_distance = j, distance
+                if (not local_unit or not osm_unit or local_unit == osm_unit) and distance < cover_distance:
+                    nearest_cover, cover_distance = j, distance
 
-            for osm_idx in nearby:
+            # Prefer marking the exact match, so e.g. a unit's own OSM node
-                osm_row = osm_gdf.iloc[osm_idx]
+            # counts as matched rather than the unit-less building around it
-                distance = local_point.distance(osm_row.geometry)
+            cover = best_match if best_match is not None else nearest_cover
+            if cover is not None:
+                matched_osm.add(cover)
 
-                local_num = self._normalize_multi(local_row.get('addr:housenumber', ''))
+            if best_match is not None:
-                osm_num   = self._normalize_multi(osm_row.get('addr:housenumber', ''))
-
-                local_street = self._normalize_street_name(str(local_row.get('addr:street', '')))
-                osm_street   = self._normalize_street_name(str(osm_row.get('addr:street', '')))
-                street_match = local_street == osm_street and local_street != ''
-
-                local_unit = self._normalize_multi(local_row.get('addr:unit'))
-                osm_unit   = self._normalize_multi(osm_row.get('addr:unit'))
-
-                if local_unit and osm_unit:
-                    unit_match = local_unit == osm_unit
-                elif local_unit or osm_unit:
-                    unit_match = False
-                else:
-                    unit_match = True
-
-                if local_num == osm_num and street_match and unit_match and distance < min_distance:
-                    min_distance = distance
-                    best_match = osm_idx
-
-            distance_m = min_distance * 111000.0
-
-            if best_match is not None and distance_m <= self.tolerance_meters:
                 existing_addresses.append({
                     'geometry': local_point,
-                    'osm_geometry': osm_gdf.iloc[best_match].geometry,
+                    'osm_geometry': osm_geoms[best_match],
-                    'local_data': dict(local_row.drop('geometry')),
+                    'local_data': local_props[i],
-                    'osm_data': dict(osm_gdf.iloc[best_match].drop('geometry')),
+                    'osm_data': osm_props[best_match],
-                    'distance_meters': distance_m,
+                    'distance_meters': min_distance * 111000.0,
                 })
             else:
-                props = dict(local_row.drop('geometry'))
+                props = dict(local_props[i])
                 props['status'] = 'new'
                 new_addresses.append({'geometry': local_point, **props})
 
-        # Second pass: find OSM addresses with no local match (potentially removed)
-        matched_osm = set()
-        for _, local_row in local_gdf.iterrows():
-            local_point = local_row.geometry
-            nearby = osm_index.query(local_point.buffer(self.tolerance_deg))
-
-            local_num    = self._normalize_multi(local_row.get('addr:housenumber', ''))
-            local_street = self._normalize_street_name(str(local_row.get('addr:street', '')))
-            local_unit   = self._normalize_multi(local_row.get('addr:unit'))
-
-            for osm_idx in nearby:
-                osm_row  = osm_gdf.iloc[osm_idx]
-                dist_m   = local_point.distance(osm_row.geometry) * 111000.0
-                osm_num  = self._normalize_multi(osm_row.get('addr:housenumber', ''))
-                osm_st   = self._normalize_street_name(str(osm_row.get('addr:street', '')))
-                osm_unit = self._normalize_multi(osm_row.get('addr:unit'))
-
-                unit_match = True
-                if local_unit and osm_unit:
-                    unit_match = local_unit == osm_unit
-
-                if (dist_m <= self.tolerance_meters and
-                        local_num == osm_num and
-                        local_street == osm_st and local_street != '' and
-                        unit_match):
-                    matched_osm.add(osm_idx)
-                    break
-
         removed_addresses = []
-        for idx, osm_row in osm_gdf.iterrows():
+        for j, osm_point in enumerate(osm_geoms):
-            if idx not in matched_osm:
+            if j not in matched_osm:
-                props = dict(osm_row.drop('geometry'))
+                props = dict(osm_props[j])
                 # geopandas reads JSON int arrays (osm_nodes) as numpy arrays,
                 # which would be written back out as a string
                 for k, v in props.items():
                     if isinstance(v, np.ndarray):
                         props[k] = [int(x) for x in v.tolist()]
                 props['status'] = 'removed'
-                removed_addresses.append({'geometry': osm_row.geometry, **props})
+                removed_addresses.append({'geometry': osm_point, **props})
 
         self._flag_close_conflicts(new_addresses)
         renames = self._flag_possible_renames(new_addresses, removed_addresses)

diff-highways.py

@@ -30,7 +30,6 @@ import warnings
 import multiprocessing as mp
 from functools import partial
 import numpy as np
-from concurrent.futures import ProcessPoolExecutor, as_completed
 import gc
 
 import importlib  
@@ -72,7 +71,6 @@ class RoadComparator:
         self.tolerance_deg = tolerance_feet / 364000.0
         self.min_gap_length_deg = min_gap_length_feet / 364000.0
         
-        print(f"Using {self.n_jobs} parallel processes with chunk size {self.chunk_size}")
         if self.exclude_unnamed:
             print("Excluding unnamed features from coverage calculation")
 
@@ -132,60 +130,27 @@ class RoadComparator:
         except Exception as e:
             raise Exception(f"Error loading {filepath}: {str(e)}")
     
-    def create_buffered_union_optimized(self, gdf: gpd.GeoDataFrame) -> Any:
+    def create_buffer_index(self, gdf: gpd.GeoDataFrame) -> Tuple[np.ndarray, STRtree]:
-        """Create a buffered union using chunked processing for memory efficiency."""
+        """Buffer every geometry and index the buffers.
-        print("Creating optimized buffered union...")
 
-        # Process in chunks to manage memory - extract geometries as lists
+        Each road is later compared only against the union of the buffers that
-        chunks = [gdf.iloc[i:i+self.chunk_size].geometry.tolist() for i in range(0, len(gdf), self.chunk_size)]
+        touch it (see _local_cover), which gives the same difference as using
-        chunk_unions = []
+        one county-wide union - buffers that don't touch a line can't change
-        
+        its difference - without every comparison paying for the size of the
-        # Use partial function for multiprocessing
+        whole county."""
-        buffer_func = partial(self._buffer_chunk, tolerance=self.tolerance_deg)
+        print("Buffering and indexing geometries...")
-        
+        buffers = np.asarray(gdf.geometry.buffer(self.tolerance_deg).values, dtype=object)
-        with ProcessPoolExecutor(max_workers=self.n_jobs) as executor:
+        return buffers, STRtree(buffers)
-            # Submit all chunk processing jobs
-            future_to_chunk = {executor.submit(buffer_func, chunk): i 
-                              for i, chunk in enumerate(chunks)}
-            
-            # Collect results as they complete
-            for future in as_completed(future_to_chunk):
-                chunk_idx = future_to_chunk[future]
-                try:
-                    chunk_union = future.result(timeout=300)  # 5 minute timeout
-                    if chunk_union and not chunk_union.is_empty:
-                        chunk_unions.append(chunk_union)
-                    print(f"Processed chunk {chunk_idx + 1}/{len(chunks)}")
-                except Exception as e:
-                    print(f"Error processing chunk {chunk_idx}: {str(e)}")
-                    # Continue processing other chunks instead of failing completely
-        
-        # Union all chunk results
-        print("Combining chunk unions...")
-        if chunk_unions:
-            final_union = unary_union(chunk_unions)
-            # Force garbage collection
-            del chunk_unions, chunks
-            gc.collect()
-            return final_union
-        else:
-            raise Exception("No valid geometries to create union")
 
     @staticmethod
-    def _buffer_chunk(geometries: List, tolerance: float) -> Any:
+    def _local_cover(geom, buffers: np.ndarray, tree: STRtree):
-        """Buffer geometries in a chunk and return their union."""
+        """Union of the target buffers that touch geom, or None if none do."""
-        try:
+        idx = tree.query(geom, predicate='intersects')
-            # Buffer all geometries in the chunk
+        if len(idx) == 0:
-            buffered = [geom.buffer(tolerance) for geom in geometries]
-            
-            # Create union of buffered geometries
-            if len(buffered) == 1:
-                return buffered[0]
-            else:
-                return unary_union(buffered)
-        except Exception as e:
-            print(f"Error in chunk processing: {str(e)}")
             return None
+        if len(idx) == 1:
+            return buffers[idx[0]]
+        return unary_union(list(buffers[idx]))
 
     def create_spatial_index(self, gdf: gpd.GeoDataFrame) -> STRtree:
         """Create spatial index for fast intersection queries."""
@@ -195,56 +160,29 @@ class RoadComparator:
         return STRtree(geometries)
     
     def find_removed_segments_optimized(self, source_gdf: gpd.GeoDataFrame,
-                                       target_union: Any) -> List[Dict[str, Any]]:
+                                       target: Tuple[np.ndarray, STRtree]) -> List[Dict[str, Any]]:
         """
-        Find segments in source_gdf that are not covered by target_union (removed roads).
+        Find segments in source_gdf that are not covered by the target buffers (removed roads).
-        Optimized with parallel processing.
         """
         print("Finding removed segments...")
+        rows = self._rows_with_cover(source_gdf, target, "removed")
+        return self._process_removed_chunk(rows, self.min_gap_length_deg)
 
-        # Split into chunks for parallel processing - convert to serializable format
+    def _rows_with_cover(self, gdf: gpd.GeoDataFrame, target: Tuple[np.ndarray, STRtree],
-        chunks = []
+                         label: str) -> List[Dict[str, Any]]:
-        for i in range(0, len(source_gdf), self.chunk_size):
+        buffers, tree = target
-            chunk_gdf = source_gdf.iloc[i:i+self.chunk_size]
+        geoms = gdf.geometry.tolist()
-            chunk_data = []
+        props = gdf.drop(columns='geometry').to_dict('records')
-            for idx, row in chunk_gdf.iterrows():
+        rows = []
-                chunk_data.append({
+        for i, (geom, p) in enumerate(zip(geoms, props)):
-                    'geometry': row.geometry,
+            rows.append({'geometry': geom, 'properties': p,
-                    'properties': dict(row.drop('geometry'))
+                         'cover': self._local_cover(geom, buffers, tree)})
-                })
+            if (i + 1) % 10000 == 0:
-            chunks.append(chunk_data)
+                print(f"  {label}: compared {i + 1}/{len(geoms)}")
-        
+        return rows
-        all_removed = []
-        
-        # Use partial function for multiprocessing
-        process_func = partial(self._process_removed_chunk, 
-                              target_union=target_union,
-                              min_length_deg=self.min_gap_length_deg)
-        
-        with ProcessPoolExecutor(max_workers=self.n_jobs) as executor:
-            # Submit all chunk processing jobs
-            future_to_chunk = {executor.submit(process_func, chunk): i 
-                              for i, chunk in enumerate(chunks)}
-            
-            # Collect results as they complete
-            for future in as_completed(future_to_chunk):
-                chunk_idx = future_to_chunk[future]
-                try:
-                    chunk_removed = future.result(timeout=300)  # 5 minute timeout
-                    all_removed.extend(chunk_removed)
-                    print(f"Processed removed chunk {chunk_idx + 1}/{len(chunks)}")
-                except Exception as e:
-                    print(f"Error processing removed chunk {chunk_idx}: {str(e)}")
-                    # Continue processing other chunks instead of failing completely
-        
-        # Clean up memory
-        del chunks
-        gc.collect()
-        return all_removed
 
     @staticmethod
-    def _process_removed_chunk(chunk_data: List[Dict], target_union: Any, 
+    def _process_removed_chunk(chunk_data: List[Dict], min_length_deg: float) -> List[Dict[str, Any]]:
-                              min_length_deg: float) -> List[Dict[str, Any]]:
         """Process a chunk of geometries to find removed segments."""
         removed_segments = []
         
@@ -261,7 +199,8 @@ class RoadComparator:
             for line in lines:
                 try:
                     # Find parts of the line that don't intersect with target_union
-                    uncovered = line.difference(target_union)
+                    cover = row_data['cover']
+                    uncovered = line.difference(cover) if cover is not None else line
                     
                     if uncovered.is_empty:
                         continue
@@ -293,56 +232,16 @@ class RoadComparator:
         return removed_segments
     
     def find_added_roads_optimized(self, source_gdf: gpd.GeoDataFrame,
-                                  target_union: Any) -> List[Dict[str, Any]]:
+                                  target: Tuple[np.ndarray, STRtree]) -> List[Dict[str, Any]]:
         """
-        Find entire roads in source_gdf that don't significantly overlap with target_union.
+        Find entire roads in source_gdf that don't significantly overlap with the target buffers.
-        Optimized with parallel processing.
         """
         print("Finding added roads...")
-        
+        rows = self._rows_with_cover(source_gdf, target, "added")
-        # Split into chunks for parallel processing - convert to serializable format
+        return self._process_added_chunk(rows, self.min_gap_length_deg)
-        chunks = []
-        for i in range(0, len(source_gdf), self.chunk_size):
-            chunk_gdf = source_gdf.iloc[i:i+self.chunk_size]
-            chunk_data = []
-            for idx, row in chunk_gdf.iterrows():
-                chunk_data.append({
-                    'geometry': row.geometry,
-                    'properties': dict(row.drop('geometry'))
-                })
-            chunks.append(chunk_data)
-        
-        all_added = []
-        
-        # Use partial function for multiprocessing
-        process_func = partial(self._process_added_chunk, 
-                              target_union=target_union,
-                              min_length_deg=self.min_gap_length_deg)
-        
-        with ProcessPoolExecutor(max_workers=self.n_jobs) as executor:
-            # Submit all chunk processing jobs
-            future_to_chunk = {executor.submit(process_func, chunk): i 
-                              for i, chunk in enumerate(chunks)}
-            
-            # Collect results as they complete
-            for future in as_completed(future_to_chunk):
-                chunk_idx = future_to_chunk[future]
-                try:
-                    chunk_added = future.result(timeout=300)  # 5 minute timeout
-                    all_added.extend(chunk_added)
-                    print(f"Processed added chunk {chunk_idx + 1}/{len(chunks)}")
-                except Exception as e:
-                    print(f"Error processing added chunk {chunk_idx}: {str(e)}")
-                    # Continue processing other chunks instead of failing completely
-        
-        # Clean up memory
-        del chunks
-        gc.collect()
-        return all_added
 
     @staticmethod
-    def _process_added_chunk(chunk_data: List[Dict], target_union: Any, 
+    def _process_added_chunk(chunk_data: List[Dict], min_length_deg: float) -> List[Dict[str, Any]]:
-                            min_length_deg: float) -> List[Dict[str, Any]]:
         """Process a chunk of geometries to find added roads."""
         added_roads = []
         
@@ -352,7 +251,8 @@ class RoadComparator:
             
             try:
                 # Check what portion of the road is not covered
-                uncovered = geom.difference(target_union)
+                cover = row_data['cover']
+                uncovered = geom.difference(cover) if cover is not None else geom
                 
                 if not uncovered.is_empty:
                     # Calculate what percentage of the original road is uncovered
@@ -465,7 +365,6 @@ class RoadComparator:
         print(f"Comparing {file1_path} and {file2_path}")
         print(f"Tolerance: {self.tolerance_feet} feet")
         print(f"Minimum significant length: {self.min_gap_length_feet} feet")
-        print(f"Parallel processing: {self.n_jobs} workers")
         print("-" * 50)
 
         # Load both files
@@ -493,20 +392,14 @@ class RoadComparator:
             print(f"Warning: CRS mismatch. Converting {file2_path} to match {file1_path}")
             gdf2 = gdf2.to_crs(gdf1.crs)
         
-        print("Creating optimized spatial unions...")
+        index1 = self.create_buffer_index(gdf1)
+        index2 = self.create_buffer_index(gdf2)
 
-        # Create buffered unions using optimized method
+        removed_roads = self.find_removed_segments_optimized(gdf1, index2)
-        union1 = self.create_buffered_union_optimized(gdf1)
+        added_roads = self.find_added_roads_optimized(gdf2, index1)
-        union2 = self.create_buffered_union_optimized(gdf2)
-        
-        print("Finding removed and added roads with parallel processing...")
-        
-        # Find roads using optimized parallel methods
-        removed_roads = self.find_removed_segments_optimized(gdf1, union2)
-        added_roads = self.find_added_roads_optimized(gdf2, union1)
 
         # Clean up memory
-        del gdf1, gdf2, union1, union2
+        del gdf1, gdf2, index1, index2
         gc.collect()
         
         return removed_roads, added_roads
@@ -635,9 +528,9 @@ Examples:
                        help='Minimum length in feet for gaps/extras to be considered significant (default: 100)')
     parser.add_argument('--output', '-o', help='Output GeoJSON file for results (optional)')
     parser.add_argument('--jobs', '-j', type=int, default=None,
-                       help='Number of parallel processes (default: CPU count - 1)')
+                       help='Ignored (kept for compatibility)')
     parser.add_argument('--chunk-size', '-c', type=int, default=1000,
-                       help='Number of geometries to process per chunk (default: 1000)')
+                       help='Ignored (kept for compatibility)')
     parser.add_argument('--exclude-unnamed', '-e', action='store_true',
                        help='Exclude features without name tags from coverage calculation (helps detect roads covered by unnamed geometry)')
 

web/server.py

@@ -3,7 +3,8 @@
 Flask web server for The Villages Import Tools
 """
 from functools import wraps
-from flask import Flask, render_template, jsonify, request, send_from_directory, session, redirect, url_for, make_response
+from flask import Flask, render_template, jsonify, request, send_from_directory, send_file, session, redirect, url_for, make_response
+from werkzeug.utils import safe_join
 from werkzeug.security import generate_password_hash, check_password_hash
 import subprocess
 import os
@@ -14,6 +15,7 @@ import shutil
 import threading
 import json
 import base64
+import gzip
 import logging
 from datetime import datetime, timezone
 import yaml
@@ -777,6 +779,9 @@ def clear_data():
                 try:
                     os.remove(fpath)
                     removed.append(f'{county["id"]}/{fname}')
+                    gz_copy = os.path.join(GZIP_CACHE_DIR, county['id'], fname + '.gz')
+                    if os.path.exists(gz_copy):
+                        os.remove(gz_copy)
                 except Exception as e:
                     errors.append(str(e))
     if errors:
@@ -820,9 +825,34 @@ def list_files():
 
     return jsonify(files)
 
+GZIP_CACHE_DIR = '/data/.gzcache'
+_gzip_lock = threading.Lock()
+
+
+def _gzipped_copy(path, filename):
+    """Path of a gzipped copy of path, (re)built when the source is newer.
+    GeoJSON compresses ~10x, which matters for 50-150 MB county files."""
+    gz_path = os.path.join(GZIP_CACHE_DIR, filename + '.gz')
+    with _gzip_lock:
+        if not os.path.exists(gz_path) or os.path.getmtime(gz_path) < os.path.getmtime(path):
+            os.makedirs(os.path.dirname(gz_path), exist_ok=True)
+            tmp = gz_path + '.tmp'
+            with open(path, 'rb') as src, gzip.open(tmp, 'wb', compresslevel=5) as dst:
+                shutil.copyfileobj(src, dst, 1024 * 1024)
+            os.replace(tmp, gz_path)
+    return gz_path
+
+
 @app.route('/data/<path:filename>')
 def serve_data(filename):
-    """Serve files from /data"""
+    """Serve files from /data (GeoJSON gzipped when the browser accepts it)"""
+    path = safe_join('/data', filename)
+    if (path and filename.endswith('.geojson') and os.path.isfile(path)
+            and 'gzip' in request.headers.get('Accept-Encoding', '')):
+        resp = make_response(send_file(_gzipped_copy(path, filename), mimetype='application/geo+json'))
+        resp.headers['Content-Encoding'] = 'gzip'
+        resp.headers['Vary'] = 'Accept-Encoding'
+        return resp
     return send_from_directory('/data', filename)
 
 # ── Config file download / upload ──────────────────────────────────────────────

web/static/map.js

@@ -123,6 +123,121 @@ function paneRenderer(pane) {
     return paneRenderers[pane];
 }
 
+// A GeoJSON layer that only creates and draws the features in (or near) the
+// current view, and nothing below `minZoom`. County and OSM files hold
+// hundreds of thousands of features, far more than a browser can draw and
+// hit-test on every pan. Layers are created on first view and then cached,
+// so a selection keeps pointing at the same layer object.
+const ViewportGeoJSON = L.GeoJSON.extend({
+    initialize: function(geojson, options) {
+        L.GeoJSON.prototype.initialize.call(this, null, options);
+        const filter = this.options.filter;
+        this._features = (geojson.features || []).filter(f => f && f.geometry && (!filter || filter(f)));
+        this._cache = new Map();   // feature -> layer
+        this._shown = new Set();   // features currently added
+        this._buildIndex();
+    },
+
+    _cellSize: 0.02, // degrees (~2 km)
+
+    _buildIndex: function() {
+        this._grid = new Map();
+        const size = this._cellSize;
+        this._features.forEach(f => {
+            let minX = Infinity, minY = Infinity, maxX = -Infinity, maxY = -Infinity;
+            const visit = c => {
+                if (typeof c[0] === 'number') {
+                    if (c[0] < minX) minX = c[0]; if (c[0] > maxX) maxX = c[0];
+                    if (c[1] < minY) minY = c[1]; if (c[1] > maxY) maxY = c[1];
+                } else c.forEach(visit);
+            };
+            visit(f.geometry.coordinates);
+            if (minX === Infinity) return;
+            for (let x = Math.floor(minX / size); x <= Math.floor(maxX / size); x++) {
+                for (let y = Math.floor(minY / size); y <= Math.floor(maxY / size); y++) {
+                    const key = x + ',' + y;
+                    let cell = this._grid.get(key);
+                    if (!cell) this._grid.set(key, cell = []);
+                    cell.push(f);
+                }
+            }
+        });
+    },
+
+    onAdd: function(map) {
+        L.GeoJSON.prototype.onAdd.call(this, map);
+        map.on('moveend', this._update, this);
+        this._update();
+    },
+
+    onRemove: function(map) {
+        map.off('moveend', this._update, this);
+        this._setHint(false);
+        L.GeoJSON.prototype.onRemove.call(this, map);
+    },
+
+    _setHint: function(show) {
+        const item = document.querySelector(`.layer-item[data-layer="${this.options.layerName}"]`);
+        if (!item) return;
+        let hint = item.querySelector('.zoom-hint');
+        if (!hint) {
+            hint = document.createElement('small');
+            hint.className = 'zoom-hint';
+            hint.style.cssText = 'margin-left:4px;color:#888;';
+            item.appendChild(hint);
+        }
+        hint.textContent = show ? '(zoom in to show)' : '';
+    },
+
+    _layerFor: function(feature) {
+        let layer = this._cache.get(feature);
+        if (layer) {
+            // Pick up style changes made while it was out of view (e.g. imported),
+            // but leave a selected layer's highlight alone.
+            if (!selectedLayers.some(s => s.layer === layer)) this.resetStyle(layer);
+            return layer;
+        }
+        // Same steps as L.GeoJSON.addData
+        layer = L.GeoJSON.geometryToLayer(feature, this.options);
+        if (!layer) return null;
+        layer.feature = L.GeoJSON.asFeature(feature);
+        layer.defaultOptions = layer.options;
+        this.resetStyle(layer);
+        if (this.options.onEachFeature) this.options.onEachFeature(feature, layer);
+        this._cache.set(feature, layer);
+        return layer;
+    },
+
+    _update: function() {
+        if (!this._map) return;
+        const tooFar = this._map.getZoom() < this.options.minZoom;
+        this._setHint(tooFar);
+        const want = new Set();
+        if (!tooFar) {
+            const b = this._map.getBounds().pad(0.25), size = this._cellSize;
+            for (let x = Math.floor(b.getWest() / size); x <= Math.floor(b.getEast() / size); x++) {
+                for (let y = Math.floor(b.getSouth() / size); y <= Math.floor(b.getNorth() / size); y++) {
+                    (this._grid.get(x + ',' + y) || []).forEach(f => want.add(f));
+                }
+            }
+        }
+        this._shown.forEach(f => {
+            if (!want.has(f)) { this.removeLayer(this._cache.get(f)); this._shown.delete(f); }
+        });
+        want.forEach(f => {
+            if (this._shown.has(f)) return;
+            const layer = this._layerFor(f);
+            if (layer) { this.addLayer(layer); this._shown.add(f); }
+        });
+    }
+});
+
+// Below these zooms the OSM/county layers aren't drawn at all.
+function viewportMinZoom(geojson) {
+    const first = (geojson.features || []).find(f => f && f.geometry);
+    return first && first.geometry.type === 'Point' ? 15 : 13;
+}
+
 // Mobile controls state
 let loadedCounty = '';
 let loadedDataType = '';
@@ -355,44 +470,21 @@ function clearSelection() {
 
 // Calculate bounds for all loaded layers
 function calculateBounds() {
-    const bounds = L.latLngBounds([]);
+    // Scan coordinates directly; building throwaway Leaflet layers for
-    let hasData = false;
+    // hundreds of thousands of features just to measure them is slow.
-
+    let minLat = Infinity, minLng = Infinity, maxLat = -Infinity, maxLng = -Infinity;
-    if (osmData && osmData.features.length > 0) {
+    const visit = c => {
-        L.geoJSON(osmData).eachLayer(layer => {
+        if (typeof c[0] === 'number') {
-            if (layer.getBounds) {
+            if (c[1] < minLat) minLat = c[1]; if (c[1] > maxLat) maxLat = c[1];
-                bounds.extend(layer.getBounds());
+            if (c[0] < minLng) minLng = c[0]; if (c[0] > maxLng) maxLng = c[0];
-            } else if (layer.getLatLng) {
+        } else c.forEach(visit);
-                bounds.extend(layer.getLatLng());
+    };
-            }
+    [osmData, diffData, countyData].forEach(data => {
-        });
+        if (!data) return;
-        hasData = true;
+        data.features.forEach(f => { if (f && f.geometry && f.geometry.coordinates) visit(f.geometry.coordinates); });
-    }
+    });
-
+    if (minLat !== Infinity) {
-    if (diffData && diffData.features.length > 0) {
+        map.fitBounds([[minLat, minLng], [maxLat, maxLng]], { padding: [50, 50] });
-        L.geoJSON(diffData).eachLayer(layer => {
-            if (layer.getBounds) {
-                bounds.extend(layer.getBounds());
-            } else if (layer.getLatLng) {
-                bounds.extend(layer.getLatLng());
-            }
-        });
-        hasData = true;
-    }
-
-    if (countyData && countyData.features.length > 0) {
-        L.geoJSON(countyData).eachLayer(layer => {
-            if (layer.getBounds) {
-                bounds.extend(layer.getBounds());
-            } else if (layer.getLatLng) {
-                bounds.extend(layer.getLatLng());
-            }
-        });
-        hasData = true;
-    }
-
-    if (hasData && bounds.isValid()) {
-        map.fitBounds(bounds, { padding: [50, 50] });
     }
 }
 
@@ -537,13 +629,15 @@ function createOsmLayer() {
 
     if (!osmData) return;
 
-    osmLayer = L.geoJSON(osmData, {
+    osmLayer = new ViewportGeoJSON(osmData, {
+        layerName: 'osm',
+        minZoom: viewportMinZoom(osmData),
         style: osmStyle,
         filter: shouldShowOsmFeature,
         pane: 'osmPane',
         renderer: paneRenderer('osmPane'),
         pointToLayer: function(feature, latlng) {
-            return L.circleMarker(latlng, osmMarkerStyle(feature));
+            return L.circleMarker(latlng, { ...osmMarkerStyle(feature), pane: 'osmPane', renderer: paneRenderer('osmPane') });
         },
         onEachFeature: function(feature, layer) {
             const isPoint = feature.geometry.type === 'Point';
@@ -642,7 +736,7 @@ function createDiffLayer() {
         pane: 'diffPane',
         renderer: paneRenderer('diffPane'),
         pointToLayer: function(feature, latlng) {
-            return L.circleMarker(latlng, diffMarkerStyle(feature));
+            return L.circleMarker(latlng, { ...diffMarkerStyle(feature), pane: 'diffPane', renderer: paneRenderer('diffPane') });
         },
         onEachFeature: function(feature, layer) {
             const isPoint = feature.geometry.type === 'Point';
@@ -713,13 +807,15 @@ function createCountyLayer() {
 
     if (!countyData) return;
 
-    countyLayer = L.geoJSON(countyData, {
+    countyLayer = new ViewportGeoJSON(countyData, {
+        layerName: 'county',
+        minZoom: viewportMinZoom(countyData),
         style: countyStyle,
         filter: shouldShowCountyFeature,
         pane: 'countyPane',
         renderer: paneRenderer('countyPane'),
         pointToLayer: function(feature, latlng) {
-            return L.circleMarker(latlng, countyMarkerStyle(feature));
+            return L.circleMarker(latlng, { ...countyMarkerStyle(feature), pane: 'countyPane', renderer: paneRenderer('countyPane') });
         },
         onEachFeature: function(feature, layer) {
             const isPoint = feature.geometry.type === 'Point';