From 4759424e853fd251b8861c9ba6c77c047c3ec337 Mon Sep 17 00:00:00 2001 From: zyphlar Date: Wed, 21 Oct 2026 19:46:02 -0700 Subject: [PATCH] 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 | 141 +++++++++++++-------------- diff-highways.py | 227 ++++++++++++------------------------------- web/server.py | 34 ++++++- web/static/map.js | 182 ++++++++++++++++++++++++++-------- 4 files changed, 296 insertions(+), 288 deletions(-) diff --git a/compare-addresses.py b/compare-addresses.py index c36802b..cda997b 100644 --- a/compare-addresses.py +++ b/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) - radius_deg = radius_meters / 111000.0 + by_num: Dict[str, List[int]] = {} + 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: - continue - osm_street = self._normalize_street_name(str(osm.get('addr:street', ''))) + if j not in street_cache: + street_cache[j] = 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(): - local_point = local_row.geometry - nearby = osm_index.query(local_point.buffer(self.tolerance_deg)) + for i, local_point in enumerate(local_geoms): + street = local_streets[i] + candidates = osm_by_key.get((local_nums[i], street), []) if street else [] + local_unit = local_units[i] - best_match = None - min_distance = float('inf') + # Existing: nearest OSM address with the same unit (or both + # 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: - osm_row = osm_gdf.iloc[osm_idx] - distance = local_point.distance(osm_row.geometry) + # Prefer marking the exact match, so e.g. a unit's own OSM node + # counts as matched rather than the unit-less building around it + 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', '')) - 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: + if best_match is not None: existing_addresses.append({ 'geometry': local_point, - 'osm_geometry': osm_gdf.iloc[best_match].geometry, - 'local_data': dict(local_row.drop('geometry')), - 'osm_data': dict(osm_gdf.iloc[best_match].drop('geometry')), - 'distance_meters': distance_m, + 'osm_geometry': osm_geoms[best_match], + 'local_data': local_props[i], + 'osm_data': osm_props[best_match], + '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(): - if idx not in matched_osm: - props = dict(osm_row.drop('geometry')) + for j, osm_point in enumerate(osm_geoms): + if j not in matched_osm: + 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 --git a/diff-highways.py b/diff-highways.py index 0dd3231..7d84cac 100644 --- a/diff-highways.py +++ b/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,61 +130,28 @@ class RoadComparator: except Exception as e: raise Exception(f"Error loading {filepath}: {str(e)}") - def create_buffered_union_optimized(self, gdf: gpd.GeoDataFrame) -> Any: - """Create a buffered union using chunked processing for memory efficiency.""" - print("Creating optimized buffered union...") - - # Process in chunks to manage memory - extract geometries as lists - chunks = [gdf.iloc[i:i+self.chunk_size].geometry.tolist() for i in range(0, len(gdf), self.chunk_size)] - chunk_unions = [] - - # Use partial function for multiprocessing - buffer_func = partial(self._buffer_chunk, tolerance=self.tolerance_deg) - - with ProcessPoolExecutor(max_workers=self.n_jobs) as executor: - # 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") - + def create_buffer_index(self, gdf: gpd.GeoDataFrame) -> Tuple[np.ndarray, STRtree]: + """Buffer every geometry and index the buffers. + + Each road is later compared only against the union of the buffers that + touch it (see _local_cover), which gives the same difference as using + 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 + whole county.""" + print("Buffering and indexing geometries...") + buffers = np.asarray(gdf.geometry.buffer(self.tolerance_deg).values, dtype=object) + return buffers, STRtree(buffers) + @staticmethod - def _buffer_chunk(geometries: List, tolerance: float) -> Any: - """Buffer geometries in a chunk and return their union.""" - try: - # Buffer all geometries in the chunk - 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)}") + def _local_cover(geom, buffers: np.ndarray, tree: STRtree): + """Union of the target buffers that touch geom, or None if none do.""" + idx = tree.query(geom, predicate='intersects') + if len(idx) == 0: 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.""" print("Creating spatial index...") @@ -194,57 +159,30 @@ class RoadComparator: geometries = gdf.geometry.tolist() return STRtree(geometries) - def find_removed_segments_optimized(self, source_gdf: gpd.GeoDataFrame, - target_union: Any) -> List[Dict[str, Any]]: + def find_removed_segments_optimized(self, source_gdf: gpd.GeoDataFrame, + target: Tuple[np.ndarray, STRtree]) -> List[Dict[str, Any]]: """ - Find segments in source_gdf that are not covered by target_union (removed roads). - Optimized with parallel processing. + Find segments in source_gdf that are not covered by the target buffers (removed roads). """ print("Finding removed segments...") - - # Split into chunks for parallel processing - convert to serializable format - 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_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 - + rows = self._rows_with_cover(source_gdf, target, "removed") + return self._process_removed_chunk(rows, self.min_gap_length_deg) + + def _rows_with_cover(self, gdf: gpd.GeoDataFrame, target: Tuple[np.ndarray, STRtree], + label: str) -> List[Dict[str, Any]]: + buffers, tree = target + geoms = gdf.geometry.tolist() + props = gdf.drop(columns='geometry').to_dict('records') + rows = [] + for i, (geom, p) in enumerate(zip(geoms, props)): + rows.append({'geometry': geom, 'properties': p, + 'cover': self._local_cover(geom, buffers, tree)}) + if (i + 1) % 10000 == 0: + print(f" {label}: compared {i + 1}/{len(geoms)}") + return rows + @staticmethod - def _process_removed_chunk(chunk_data: List[Dict], target_union: Any, - min_length_deg: float) -> List[Dict[str, Any]]: + def _process_removed_chunk(chunk_data: List[Dict], 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 @@ -292,57 +231,17 @@ class RoadComparator: return removed_segments - def find_added_roads_optimized(self, source_gdf: gpd.GeoDataFrame, - target_union: Any) -> List[Dict[str, Any]]: + def find_added_roads_optimized(self, source_gdf: gpd.GeoDataFrame, + target: Tuple[np.ndarray, STRtree]) -> List[Dict[str, Any]]: """ - Find entire roads in source_gdf that don't significantly overlap with target_union. - Optimized with parallel processing. + Find entire roads in source_gdf that don't significantly overlap with the target buffers. """ print("Finding added roads...") - - # Split into chunks for parallel processing - convert to serializable format - 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 - + rows = self._rows_with_cover(source_gdf, target, "added") + return self._process_added_chunk(rows, self.min_gap_length_deg) + @staticmethod - def _process_added_chunk(chunk_data: List[Dict], target_union: Any, - min_length_deg: float) -> List[Dict[str, Any]]: + def _process_added_chunk(chunk_data: List[Dict], 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...") - - # Create buffered unions using optimized method - union1 = self.create_buffered_union_optimized(gdf1) - 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) - + index1 = self.create_buffer_index(gdf1) + index2 = self.create_buffer_index(gdf2) + + removed_roads = self.find_removed_segments_optimized(gdf1, index2) + added_roads = self.find_added_roads_optimized(gdf2, index1) + # 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)') diff --git a/web/server.py b/web/server.py index a2d0dac..8466f36 100644 --- a/web/server.py +++ b/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/') 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 ────────────────────────────────────────────── diff --git a/web/static/map.js b/web/static/map.js index 534d8c6..112754c 100644 --- a/web/static/map.js +++ b/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([]); - let hasData = false; - - if (osmData && osmData.features.length > 0) { - L.geoJSON(osmData).eachLayer(layer => { - if (layer.getBounds) { - bounds.extend(layer.getBounds()); - } else if (layer.getLatLng) { - bounds.extend(layer.getLatLng()); - } - }); - hasData = true; - } - - if (diffData && diffData.features.length > 0) { - 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] }); + // Scan coordinates directly; building throwaway Leaflet layers for + // hundreds of thousands of features just to measure them is slow. + let minLat = Infinity, minLng = Infinity, maxLat = -Infinity, maxLng = -Infinity; + const visit = c => { + if (typeof c[0] === 'number') { + if (c[1] < minLat) minLat = c[1]; if (c[1] > maxLat) maxLat = c[1]; + if (c[0] < minLng) minLng = c[0]; if (c[0] > maxLng) maxLng = c[0]; + } else c.forEach(visit); + }; + [osmData, diffData, countyData].forEach(data => { + if (!data) return; + data.features.forEach(f => { if (f && f.geometry && f.geometry.coordinates) visit(f.geometry.coordinates); }); + }); + if (minLat !== Infinity) { + map.fitBounds([[minLat, minLng], [maxLat, maxLng]], { 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';