import numpy as np
from rtree import index
from concurrent.futures import ThreadPoolExecutorclass GeoFenceEngine:def __init__(self, fences):"""初始化地理围栏引擎"""self.idx = index.Index()self.fences = []for i, fence in enumerate(fences):# 为每个围栏创建空间索引self.idx.insert(i, fence['bbox'])self.fences.append({'poly': np.array(fence['poly']),'holes': [np.array(h) for h in fence.get('holes', [])],'attrs': fence.get('attrs', {})})def query(self, points, batch_size=1000):"""批量查询点所在的围栏"""results = []with ThreadPoolExecutor() as executor:for i in range(0, len(points), batch_size):batch = points[i:i+batch_size]results.extend(executor.map(self._query_point, batch))return resultsdef _query_point(self, point):"""查询单个点所在的围栏"""x, y = pointfor fid in self.idx.intersection((x, y, x, y)):fence = self.fences[fid]if self._point_in_polygon(point, fence['poly']):if not any(self._point_in_polygon(point, h) for h in fence['holes']):return {'fence_id': fid, 'attrs': fence['attrs']}return None@staticmethoddef _point_in_polygon(point, polygon):"""优化的点在多边形内判断"""# 实现同前文...pass# 生产环境使用示例
if __name__ == "__main__":# 定义多个复杂围栏fences = [{'poly': [(0,0),(5,0),(5,5),(0,5)],'holes': [[(1,1),(2,4),(3,1)]],'bbox': (0,0,5,5),'attrs': {'name': '中央公园', 'type': 'park'}},# 可添加更多围栏...]engine = GeoFenceEngine(fences)# 批量查询(模拟10000个点)np.random.seed(42)points = np.random.uniform(-1, 6, size=(10000, 2))import timestart = time.time()results = engine.query(points)print(f"处理{len(points)}个点耗时: {time.time()-start:.3f}秒")print("前5个结果:", results[:5])

from sklearn.cluster import OPTICS
from datetime import datetime
import matplotlib.pyplot as pltclass TrajectoryCluster:def __init__(self, spatial_eps=100, temporal_eps=300, min_samples=5):"""spatial_eps: 空间阈值(米)temporal_eps: 时间阈值(秒)min_samples: 最小聚类点数"""self.spatial_eps = spatial_eps / 6371000  # 转换为弧度self.temporal_eps = temporal_epsself.min_samples = min_samplesdef fit_predict(self, trajectories):"""执行轨迹聚类"""# 准备时空特征矩阵features = []for traj in trajectories:for point in traj['points']:# 空间坐标(转换为弧度)lat_rad = np.radians(point['lat'])lon_rad = np.radians(point['lon'])# 时间戳(转换为秒)ts = point['time'].timestamp()features.append([lat_rad, lon_rad, ts])features = np.array(features)# 自定义距离度量def spatiotemporal_dist(a, b):spatial_dist = haversine(a[:2], b[:2])time_dist = abs(a[2]-b[2])# 标准化后加权return (spatial_dist/self.spatial_eps + time_dist/self.temporal_eps)/2# 使用OPTICS算法(比DBSCAN更适合时空数据)clustering = OPTICS(min_samples=self.min_samples,metric=spatiotemporal_dist,cluster_method='xi').fit(features)return clustering.labels_def visualize(self, trajectories, labels):"""可视化聚类结果"""plt.figure(figsize=(12, 8))unique_labels = set(labels)colors = plt.cm.tab20(np.linspace(0, 1, len(unique_labels)))for label, color in zip(unique_labels, colors):if label == -1:color = 'gray'  # 噪声点mask = (labels == label)cluster_points = []for i, traj in enumerate(trajectories):if mask[i]:for point in traj['points']:plt.scatter(point['lon'], point['lat'],color=color, alpha=0.5, s=10)plt.xlabel('经度')plt.ylabel('纬度')plt.title('轨迹聚类结果')plt.grid(True)plt.show()def haversine(a, b):"""计算球面距离"""# 实现略...pass# 示例使用
if __name__ == "__main__":# 模拟轨迹数据def random_trajectory(base_lat, base_lon, base_time, n_points=10):points = []time = base_timefor _ in range(n_points):points.append({'lat': base_lat + np.random.uniform(-0.01, 0.01),'lon': base_lon + np.random.uniform(-0.01, 0.01),'time': time})time += timedelta(minutes=5)return {'points': points}from datetime import datetime, timedeltabase_time = datetime.now()trajectories = [random_trajectory(39.9, 116.4, base_time),  # 轨迹组1random_trajectory(39.9, 116.4, base_time + timedelta(hours=1)),random_trajectory(39.8, 116.3, base_time),  # 轨迹组2random_trajectory(39.8, 116.3, base_time + timedelta(hours=2)),random_trajectory(40.0, 116.5, base_time)   # 噪声轨迹]clusterer = TrajectoryCluster(spatial_eps=500, temporal_eps=1800)labels = clusterer.fit_predict(trajectories)clusterer.visualize(trajectories, labels)

import networkx as nx
from heapq import heappush, heappopclass RealTimeRouter:def __init__(self, graph):"""初始化路由引擎"""self.base_graph = graphself.current_graph = graph.copy()self.ch_graph = self._preprocess_contraction_hierarchies()def _preprocess_contraction_hierarchies(self):"""预处理收缩层次结构"""# 简化的CH预处理(实际工程需要更复杂实现)G = self.base_graph.copy()levels = {node: 0 for node in G.nodes()}# 这里应实现完整的CH预处理流程# 简化为只添加少量捷径for node in sorted(G.nodes(), key=lambda x: -G.degree(x)):neighbors = list(G.neighbors(node))for i, u in enumerate(neighbors):for v in neighbors[i+1:]:if not G.has_edge(u, v):G.add_edge(u, v, weight=G[u][node]['weight']+G[node][v]['weight'])levels[node] += 1return Gdef update_traffic(self, edge_updates):"""更新实时交通状况"""for (u, v), attrs in edge_updates.items():if self.current_graph.has_edge(u, v):self.current_graph[u][v].update(attrs)def shortest_path(self, source, target, criteria='time'):"""多标准最短路径查询"""def cost(u, v, data):# 多标准代价计算if criteria == 'time':return data.get('time', data['weight'])elif criteria == 'distance':return data.get('distance', data['weight'])elif criteria == 'composite':return 0.7*data.get('time',1) + 0.3*data.get('toll',0)return data['weight']# 使用CH加速查询return self._bidirectional_dijkstra(source, target, cost)def _bidirectional_dijkstra(self, source, target, cost_func):"""双向Dijkstra算法实现"""# 实现略...pass# 生产环境示例
if __name__ == "__main__":# 创建测试路网G = nx.DiGraph()# 添加节点(带坐标)nodes = {1: (116.404, 39.915), 2: (116.410, 39.920),3: (116.405, 39.925), 4: (116.415, 39.918)}for n, pos in nodes.items():G.add_node(n, pos=pos)# 添加边(带多种权重)edges = [(1,2, {'weight':1.5, 'time':10, 'distance':1500, 'toll':0}),(2,3, {'weight':1.2, 'time':8, 'distance':1200, 'toll':5}),(3,4, {'weight':2.0, 'time':15, 'distance':2000, 'toll':0}),(1,4, {'weight':3.0, 'time':25, 'distance':3000, 'toll':10})]G.add_edges_from(edges)router = RealTimeRouter(G)# 模拟实时交通更新router.update_traffic({(1,2): {'time':15},  # 拥堵导致时间增加(3,4): {'time':10}   # 路况改善})# 查询路径print("最快路径:", router.shortest_path(1, 4, 'time'))print("最短距离:", router.shortest_path(1, 4, 'distance'))print("综合最优:", router.shortest_path(1, 4, 'composite'))

import dask.dataframe as dd
from dask.distributed import Client
from quadtree import QuadTreeclass SpatialJoiner:def __init__(self, left_df, right_df, spatial_cols=['lon','lat']):"""初始化空间连接器"""self.left = left_dfself.right = right_dfself.spatial_cols = spatial_colsself.client = Client()  # 启动Dask集群# 构建空间索引self._build_index()def _build_index(self):"""构建QuadTree索引"""# 计算空间范围all_points = dd.concat([self.left[self.spatial_cols],self.right[self.spatial_cols]]).compute()x_min, y_min = all_points.min()x_max, y_max = all_points.max()self.qt = QuadTree((x_min, y_min, x_max, y_max),max_items=1000)# 为右表构建索引right_points = self.right[self.spatial_cols].compute()for idx, row in right_points.iterrows():self.qt.insert(row.values, idx)def join(self, distance, how='inner'):"""执行空间连接"""def find_matches(partition):matches = []for idx, row in partition.iterrows():point = row[self.spatial_cols].valuesneighbors = self.qt.query_radius(point, distance)for n in neighbors:matches.append((idx, n))return pd.DataFrame(matches, columns=['left_idx','right_idx'])# 并行处理左表分区matches = self.left.map_partitions(find_matches,meta={'left_idx': 'int64', 'right_idx': 'int64'}).compute()# 执行连接left_df = self.left.compute()right_df = self.right.compute()result = pd.merge(left_df.reset_index().rename(columns={'index':'left_idx'}),matches,on='left_idx').merge(right_df.reset_index().rename(columns={'index':'right_idx'}),on='right_idx')return result.drop(columns=['left_idx','right_idx'])# 使用示例
if __name__ == "__main__":# 创建测试数据集(实际应从文件读取)n_points = 100000left_data = {'id': range(n_points),'lon': np.random.uniform(116.3, 116.5, n_points),'lat': np.random.uniform(39.8, 40.0, n_points),'value1': np.random.rand(n_points)}right_data = {'id': range(n_points//2),'lon': np.random.uniform(116.35, 116.45, n_points//2),'lat': np.random.uniform(39.85, 39.95, n_points//2),'value2': np.random.rand(n_points//2)}# 创建Dask DataFrameleft_df = dd.from_pandas(pd.DataFrame(left_data), npartitions=4)right_df = dd.from_pandas(pd.DataFrame(right_data), npartitions=2)# 执行空间连接joiner = SpatialJoiner(left_df, right_df)result = joiner.join(distance=0.01)  # 约1公里print(f"找到{len(result)}对空间邻近记录")print(result.head())

import rasterio
from rasterio.windows import Window
from rasterio import Affine
import numpy as npclass RasterProcessor:def __init__(self, file_path):"""初始化栅格处理器"""self.src = rasterio.open(file_path)self.profile = self.src.profiledef process_by_blocks(self, block_size=1024, process_func=None):"""分块处理栅格数据"""if process_func is None:process_func = lambda x: x  # 默认不处理# 准备输出文件out_profile = self.profile.copy()out_path = 'processed_' + self.src.namewith rasterio.open(out_path, 'w', **out_profile) as dst:for ji, window in self.src.block_windows():# 读取当前块block = self.src.read(window=window)# 处理数据processed = process_func(block)# 写入结果dst.write(processed, window=window)return out_pathdef zonal_stats(self, polygons, stats=['mean','max','min']):"""区域统计"""results = []for poly in polygons:# 获取多边形边界框minx, miny, maxx, maxy = poly.bounds# 转换到栅格坐标window = self.src.window(minx, miny, maxx, maxy)# 读取数据data = self.src.read(window=window)# 创建多边形掩膜from rasterio.features import geometry_maskmask = geometry_mask([poly],transform=self.src.window_transform(window),out_shape=data.shape[1:])# 应用统计masked_data = np.ma.masked_array(data, mask=mask)stat_results = {}if 'mean' in stats:stat_results['mean'] = masked_data.mean()if 'max' in stats:stat_results['max'] = masked_data.max()if 'min' in stats:stat_results['min'] = masked_data.min()results.append({'polygon_id': poly.get('id', 'unknown'),**stat_results})return results# 使用示例
if __name__ == "__main__":# 示例栅格处理def ndvi_calculator(block):"""计算NDVI植被指数"""red = block[0].astype('float32')nir = block[1].astype('float32')ndvi = (nir - red) / (nir + red + 1e-6)return ndvi[np.newaxis, ...]  # 保持三维形状processor = RasterProcessor('example.tif')processed_path = processor.process_by_blocks(process_func=ndvi_calculator)print(f"处理结果已保存到: {processed_path}")# 示例区域统计from shapely.geometry import shapepolygons = [{'id': 'field1', 'geometry': shape({'type': 'Polygon','coordinates': [[[116.4,39.9],[116.41,39.9],[116.41,39.91],[116.4,39.91]]})}]stats = processor.zonal_stats(polygons)print("区域统计结果:", stats)