diff --git a/example/Example2-facade shadow.ipynb b/example/Example2-facade shadow.ipynb new file mode 100644 index 0000000..81dba78 --- /dev/null +++ b/example/Example2-facade shadow.ipynb @@ -0,0 +1,119 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "import pandas as pd\n", + "import warnings\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "\n", + "import pybdshadow\n", + "MAPBOX_ACCESS_TOKEN = \"pk.eyJ1IjoibmkxbzEiLCJhIjoiY2t3ZDgzMmR5NDF4czJ1cm84Z3NqOGt3OSJ9.yOYP6pxDzXzhbHfyk3uORg\"\n", + "bounds = [139.803137,35.690984,139.804437,35.692684]\n", + "buildings_gdf = pybdshadow.get_buildings_by_bounds(139.804337,35.692584,139.804437,35.692684,MAPBOX_ACCESS_TOKEN)\n", + "\n", + "\n", + "buildings_gdf['x'] = buildings_gdf.centroid.x\n", + "buildings_gdf['y'] = buildings_gdf.centroid.y\n", + "buildings_gdf = buildings_gdf[(buildings_gdf['x'] > bounds[0]) &\n", + " (buildings_gdf['x'] < bounds[2]) &\n", + " (buildings_gdf['y'] > bounds[1]) &\n", + " (buildings_gdf['y'] < bounds[3])]\n", + "buildings_gdf.plot()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "precision = 3600\n", + "date = '2022-01-01'\n", + "wallsunshine = pybdshadow.cal_sunshine_facade(buildings_gdf, date, precision)\n", + "\n", + "\n", + "floorsunshine = pybdshadow.cal_sunshine(buildings_gdf,\n", + " day=date,\n", + " roof=False,\n", + " accuracy='vector',\n", + " precision=precision)\n", + "floorsunshine['height'] = 0\n", + "\n", + "roofsunshine = pybdshadow.cal_sunshine(buildings_gdf,\n", + " day=date,\n", + " roof=True,\n", + " accuracy='vector',\n", + " precision=precision\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "User Guide: https://docs.kepler.gl/docs/keplergl-jupyter\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "90b88e3da79444708001c3ebdbf20e5d", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "KeplerGl(config={'version': 'v1', 'config': {'visState': {'filters': [], 'layers': [{'id': 'lz48o4', 'type': '…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\n", + "# 生成立体建筑物\n", + "plane_sunshine = pd.concat([floorsunshine,roofsunshine],axis=0)\n", + "\n", + "plane_sunshine['geometry'] = plane_sunshine.apply(lambda x: pybdshadow.extrude_poly(x['geometry'],x['height']),axis=1)\n", + "\n", + "final_shadows_sunshinetime = pd.concat([plane_sunshine,wallsunshine],axis=0)\n", + "\n", + "vis = pybdshadow.show_sunshine(sunshine = final_shadows_sunshinetime,\n", + " zoom='auto',vis_height = 1000)\n", + "vis" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py38_native", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/src/pybdshadow/__init__.py b/src/pybdshadow/__init__.py index e2fa89d..96e1710 100644 --- a/src/pybdshadow/__init__.py +++ b/src/pybdshadow/__init__.py @@ -32,7 +32,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """ -__version__ = '0.3.4' +__version__ = '0.3.5' __author__ = 'Qing Yu ' # module level doc-string @@ -40,6 +40,10 @@ `pybdshadow` - Python package to generate building shadow geometry. """ from .pybdshadow import * +from .get_buildings import ( + get_buildings_by_polygon, + get_buildings_by_bounds, +) from .pybdshadow import ( bdshadow_sunlight, bdshadow_pointlight @@ -49,6 +53,7 @@ ) from .visualization import ( show_bdshadow, + show_sunshine, ) from .analysis import ( cal_sunshine, @@ -57,6 +62,14 @@ get_timetable ) +from .facade import ( + cal_sunshine_facade +) + +from .utils import ( + extrude_poly +) + __all__ = ['bdshadow_sunlight', 'bdshadow_pointlight', 'bd_preprocess', @@ -64,5 +77,10 @@ 'cal_sunshine', 'cal_sunshadows', 'cal_shadowcoverage', - 'get_timetable' + 'get_timetable', + 'get_buildings_by_polygon', + 'get_buildings_by_bounds', + 'cal_sunshine_facade', + 'show_sunshine', + 'extrude_poly' ] diff --git a/src/pybdshadow/analysis.py b/src/pybdshadow/analysis.py index dcfbb00..c60cfe8 100644 --- a/src/pybdshadow/analysis.py +++ b/src/pybdshadow/analysis.py @@ -1,13 +1,13 @@ import pandas as pd from suncalc import get_times -from shapely.geometry import MultiPolygon +from shapely.geometry import MultiPolygon,Polygon import transbigdata as tbd import geopandas as gpd from .pybdshadow import ( bdshadow_sunlight, ) from .preprocess import bd_preprocess - +from .utils import count_overlapping_features def get_timetable(lon, lat, dates=['2022-01-01'], precision=3600, padding=1800): # generate timetable with given interval @@ -56,6 +56,8 @@ def cal_sunshine(buildings, day='2022-01-01', roof=False, grids=gpd.GeoDataFrame grids generated by TransBigData in study area, each grids have a `time` column store the sunshine time ''' + + # calculate day time duration lon, lat = buildings['geometry'].iloc[0].bounds[:2] date = pd.to_datetime(day+' 12:45:33.959797119') @@ -73,22 +75,47 @@ def cal_sunshine(buildings, day='2022-01-01', roof=False, grids=gpd.GeoDataFrame if accuracy == 'vector': if roof: shadows = shadows[shadows['type'] == 'roof'] - shadows = bd_preprocess(shadows) - shadows = shadows.groupby(['date', 'type'])['geometry'].apply( - lambda df: MultiPolygon(list(df)).buffer(0)).reset_index() - shadows = bd_preprocess(shadows) - shadows = count_overlapping_features(shadows) + if len(shadows)>0: + shadows = bd_preprocess(shadows) + shadows = shadows.groupby(['date', 'type','height'])['geometry'].apply( + lambda df: MultiPolygon(list(df)).buffer(0)).reset_index() + shadows = bd_preprocess(shadows) + + # 额外：增加屋顶面 + shadows = pd.concat([shadows, buildings]) + #return shadows + shadows = shadows.groupby('height').apply(count_overlapping_features).reset_index() + shadows['count'] -= 1 else: shadows = shadows[shadows['type'] == 'ground'] + shadows = bd_preprocess(shadows) shadows = shadows.groupby(['date', 'type'])['geometry'].apply( - lambda df: MultiPolygon(list(df)).buffer(0)).reset_index() + lambda df: MultiPolygon(list(df)).buffer(0)).reset_index() shadows = bd_preprocess(shadows) - shadows = count_overlapping_features(shadows) + + # 额外：增加地面面 + minpos = shadows.bounds[['minx','miny']].min() + maxpos = shadows.bounds[['maxx','maxy']].max() + + ground = gpd.GeoDataFrame(geometry=[ + Polygon([ + [minpos['minx'],minpos['miny']], + [minpos['minx'],maxpos['maxy']], + [maxpos['maxx'],maxpos['maxy']], + [maxpos['maxx'],minpos['miny']], + ]) + ]) + shadows = pd.concat([shadows, + ground + ]) + shadows = count_overlapping_features(shadows,buffer=False) + shadows['count'] -= 1 + shadows['time'] = shadows['count']*precision shadows['Hour'] = sunlighthour-shadows['time']/3600 - shadows.loc[shadows['Hour'] <= 0, 'Hour'] = 0 + #shadows.loc[shadows['Hour'] <= 0, 'Hour'] = 0 return shadows else: # Grid analysis of shadow cover duration(ground). @@ -97,7 +124,7 @@ def cal_sunshine(buildings, day='2022-01-01', roof=False, grids=gpd.GeoDataFrame grids['Hour'] = sunlighthour-grids['time']/3600 return grids - + def cal_sunshadows(buildings, cityname='somecity', dates=['2022-01-01'], precision=3600, padding=1800, roof=True, include_building=True, save_shadows=False, printlog=False): @@ -222,17 +249,3 @@ def cal_shadowcoverage(shadows_input, buildings, grids=gpd.GeoDataFrame(), roof= return grids - -def count_overlapping_features(gdf): - import shapely - bounds = gdf.geometry.exterior.unary_union - new_polys = list(shapely.ops.polygonize(bounds)) - new_gdf = gpd.GeoDataFrame(geometry=new_polys) - new_gdf['id'] = range(len(new_gdf)) - new_gdf_centroid = new_gdf.copy() - new_gdf_centroid['geometry'] = new_gdf.centroid - overlapcount = gpd.sjoin(new_gdf_centroid, gdf) - overlapcount = overlapcount.groupby( - ['id'])['index_right'].count().rename('count').reset_index() - out_gdf = pd.merge(new_gdf, overlapcount) - return out_gdf diff --git a/src/pybdshadow/facade.py b/src/pybdshadow/facade.py new file mode 100644 index 0000000..0b49c17 --- /dev/null +++ b/src/pybdshadow/facade.py @@ -0,0 +1,629 @@ +from .utils import ( + lonlat2aeqd, + aeqd2lonlat_3d, + has_normal, + count_overlapping_features, + calculate_normal, + make_clockwise +) +from .analysis import get_timetable +from .pybdshadow import bdshadow_sunlight +from shapely.geometry import Polygon, MultiPolygon, MultiPolygon, GeometryCollection +import geopandas as gpd +import numpy as np +import pandas as pd +import suncalc +from suncalc import get_times + +def cal_multiple_wall_overlap_count(walls): + def to_3d(result_wall): + # 2d转3d + # 提取并集或交集的坐标 + if isinstance(result_wall, Polygon): + # 如果结果是单一多边形 + coords = np.array(result_wall.exterior.coords) + elif isinstance(result_wall, MultiPolygon): + # 如果结果是多边形集合，合并它们的坐标 + coords = np.concatenate([np.array(poly.exterior.coords) + for poly in result_wall.geoms]) + else: + return np.array([]) + # 计算平面方程中的常数 d + d = -np.dot(normal_vector, p1) + + # 根据选择的坐标轴，反推缺失的坐标 + if coords_plane == (1, 2): # y和z轴，需要反推x坐标 + if nx == 0: + raise ValueError("平面方程无法解出唯一的x值") # 无法处理垂直于x轴的平面 + x = -(ny * coords[:, 0] + nz * coords[:, 1] + d) / nx + result = np.c_[x, coords] + + elif coords_plane == (0, 2): # x和z轴，需要反推y坐标 + if ny == 0: + raise ValueError("平面方程无法解出唯一的y值") + y = -(nx * coords[:, 0] + nz * coords[:, 1] + d) / ny + result = np.c_[coords[:, 0], y, coords[:, 1]] + return Polygon(result) + + walls = list(walls['geometry'].apply(lambda x: list(x.exterior.coords))) + p1 = walls[0][0] + normal_vector = calculate_normal(walls[0]) + nx, ny, nz = normal_vector + + # 根据法向量的方向选择坐标轴 + if abs(nx) > abs(ny): + # 法向量主要指向x轴，选择y和z轴 + coords_plane = (1, 2) + else: + # 法向量主要指向y轴，选择x和z轴 + coords_plane = (0, 2) + + gdf = gpd.GeoDataFrame( + geometry=[Polygon(np.array(wall)[:, coords_plane]) for wall in walls]) + + overlap = count_overlapping_features(gdf, buffer=False) + + overlap['geometry'] = overlap['geometry'].apply(to_3d) + return overlap + +def cal_multiple_wall_union(walls): + """ + 计算多个墙（平面）的并集。 + + 此函数通过接收一个包含多个墙面顶点的列表来计算它们的并集。每个墙面由至少三个顶点在三维空间中定义。 + + 参数: + walls: 一个三维数组，其中每个元素是一个墙面的顶点列表。每个墙面是由三个或更多的三维点（x, y, z）组成的列表。 + 例如 walls = [wall1,wall2] 其中，wall1: 第一个墙的坐标点列表，格式为 [[x1, y1, z1], [x2, y2, z2], ...]，wall2: 第二个墙的坐标点列表，格式为 [[x1, y1, z1], [x2, y2, z2], ...] + + 返回: + result: 并集的坐标点数组 + + """ + + walls = np.array(walls) + p1, p2, p3 = walls[0][:3] + normal_vector = calculate_normal(walls[0]) + nx, ny, nz = normal_vector + + # 根据法向量的方向选择坐标轴 + if abs(nx) > abs(ny): + # 法向量主要指向x轴，选择y和z轴 + coords_plane = (1, 2) + else: + # 法向量主要指向y轴，选择x和z轴 + coords_plane = (0, 2) + + result_wall = MultiPolygon([Polygon(wall) + for wall in walls[:, :, coords_plane]]).buffer(0) + + # 提取并集或交集的坐标 + if isinstance(result_wall, Polygon): + # 如果结果是单一多边形 + coords = np.array(result_wall.exterior.coords) + elif isinstance(result_wall, MultiPolygon): + # 如果结果是多边形集合，合并它们的坐标 + coords = np.concatenate([np.array(poly.exterior.coords) + for poly in result_wall.geoms]) + else: + return np.array([]) + # 计算平面方程中的常数 d + d = -np.dot(normal_vector, p1) + + # 根据选择的坐标轴，反推缺失的坐标 + if coords_plane == (1, 2): # y和z轴，需要反推x坐标 + if nx == 0: + raise ValueError("平面方程无法解出唯一的x值") # 无法处理垂直于x轴的平面 + x = -(ny * coords[:, 0] + nz * coords[:, 1] + d) / nx + result = np.c_[x, coords] + + elif coords_plane == (0, 2): # x和z轴，需要反推y坐标 + if ny == 0: + raise ValueError("平面方程无法解出唯一的y值") + y = -(nx * coords[:, 0] + nz * coords[:, 1] + d) / ny + result = np.c_[coords[:, 0], y, coords[:, 1]] + + return result + +def cal_wall_overlap(wall1, wall2, method='intersection'): + """ + 计算两个墙（平面）的并集或交集。 + + 参数: + wall1: 第一个墙的坐标点列表，格式为 [[x1, y1, z1], [x2, y2, z2], ...] + wall2: 第二个墙的坐标点列表，格式为 [[x1, y1, z1], [x2, y2, z2], ...] + method: 计算方式，'union' 为并集，'intersection' 为交集，默认为 'union' + + 返回: + result: 并集或交集的坐标点数组，格式与输入格式相同 + + 思路： + 核心思路是根据两个墙面（平面）的法向量来确定它们主要垂直于哪个坐标轴。然后，选择两个与主轴垂直的坐标轴来表征这些平面。 + 例如，如果一个平面主要垂直于 x 轴，那么我们使用 y 和 z 轴。 + 基于这些坐标轴，使用 Shapely 库创建多边形来代表每个墙面。 + 接下来，计算这两个多边形在xz或yz平面的并集或交集。 + 最后，根据先前选择的坐标轴，使用平面方程反推缺失的坐标值，以获取完整的三维坐标点集。 + 这样就能够处理那些垂直于不同坐标轴的墙面，并计算它们的重叠部分。 + """ + + # 计算法向量 + plane = wall1 + p1, p2, p3 = plane[:3] + + normal_vector = calculate_normal(plane) + nx, ny, nz = normal_vector + + # 根据法向量的方向选择坐标轴 + if abs(nx) > abs(ny): + # 法向量主要指向x轴，选择y和z轴 + coords_plane = (1, 2) + else: + # 法向量主要指向y轴，选择x和z轴 + coords_plane = (0, 2) + + # 将墙的坐标转换为多边形 + poly1 = Polygon(np.array(wall1)[:, coords_plane]) + if not poly1.is_valid: + poly1 = poly1.buffer(0) + poly2 = Polygon(np.array(wall2)[:, coords_plane]) + if not poly2.is_valid: + poly2 = poly2.buffer(0) + + # 根据方法计算并集或交集 + if method == 'union': + result_wall = poly1.union(poly2) + elif method == 'intersection': + result_wall = poly1.intersection(poly2) + + # 提取并集或交集的坐标 + final_result = None + if isinstance(result_wall, Polygon): + final_result = result_wall + elif isinstance(result_wall, GeometryCollection): + # 如果是几何集合，处理每个元素 + for geom in result_wall.geoms: + if isinstance(geom, Polygon): + if final_result is None: + final_result = geom + else: + if method == 'union': + final_result = final_result.union(geom) + elif method == 'intersection': + final_result = final_result.intersection(geom) + + if final_result is not None and isinstance(final_result, Polygon): + # 提取最终结果的坐标 + coords = np.array(final_result.exterior.coords) + else: + # 如果没有有效的 Polygon，返回空数组 + return np.array([]) + # 计算平面方程中的常数 d + d = -np.dot(normal_vector, p1) + + # 根据选择的坐标轴，反推缺失的坐标 + if coords_plane == (1, 2): # y和z轴，需要反推x坐标 + if nx == 0: + raise ValueError("平面方程无法解出唯一的x值") # 无法处理垂直于x轴的平面 + x = -(ny * coords[:, 0] + nz * coords[:, 1] + d) / nx + result = np.c_[x, coords] + + elif coords_plane == (0, 2): # x和z轴，需要反推y坐标 + if ny == 0: + raise ValueError("平面方程无法解出唯一的x值", str(wall1), str(wall2)) + # raise ValueError("平面方程无法解出唯一的y值") # 无法处理垂直于y轴的平面 + return None + + y = -(nx * coords[:, 0] + nz * coords[:, 1] + d) / ny + result = np.c_[coords[:, 0], y, coords[:, 1]] + + return result + +def calculate_wall_normal_vector(wall): + """ + 计算给定墙面的法向量。 + + 假设墙面顶点按顺时针方向给出。 + + 参数: + wall : list 或 ndarray + 墙面上的点的集合，其中每个点是一个三维坐标 [x, y, z]。 + + 返回: + normal_vector : ndarray + 墙面的法向量。 + """ + # 取墙面的前三个点 + p1, p2, p3 = np.array(wall[:3]) + + # 计算两个向量，这两个向量在墙面上并且相互垂直 + v1 = p2 - p1 + v2 = p3 - p1 + + # 计算法向量，使用叉积 + normal_vector = np.cross(v1, v2) + + # 标准化法向量 + normal_vector = normal_vector / np.linalg.norm(normal_vector) + + return normal_vector + +def calculate_wall_plane(wall): + p1, p2, p3 = np.array(wall[:3]) + v1 = p2 - p1 + v2 = p3 - p1 + normal_vector = np.cross(v1, v2) + + # 计算平面方程 Ax + By + Cz + D = 0 中的 A, B, C, D + A, B, C = normal_vector + D = -np.dot(normal_vector, p1) + + return A, B, C, D + +def projection_on_wall_single(plane, sun_vec, wall, wall_normal): + """ + 计算单个墙在单个平面上由太阳光照射产生的投影点，并根据夹角判断是否保留投影点。 + + 输入: + plane : ndarray + 平面方程的系数[A, B, C, D]。 + sun_vec : list 或 ndarray + 太阳光的方向向量，格式为 [vx, vy, vz]。 + wall : ndarray + 墙面上的点的集合，格式为 4*3。 + wall_normal : list 或 ndarray + 被投影墙面的法向量。 + + 输出: + intersections : ndarray + 4*3大小的矩阵，平面每个点上由太阳光照射产生的投影点的坐标，不合理的投影点标记为None。 + """ + sun_vec = np.array(sun_vec) + wall = np.array(wall) + A, B, C, D = plane + + denominator = np.dot(np.array([A, B, C]), sun_vec) + t = -(np.dot(wall, np.array([A, B, C])) + D) / denominator + intersections = wall + np.outer(t, sun_vec) + + # 处理特殊情况，平行或无交点 + intersections[denominator == 0] = np.array([None, None, None]) + + # 检查每个投影点是否在墙面法向量的“背面” + for i in range(len(intersections)): + + if not np.any(np.isnan(intersections[i])): + + # 计算从墙面顶点到投影点的向量 + vector_to_projection = intersections[i] - wall[i] + if np.linalg.norm(vector_to_projection) > 0: + + # 计算夹角 + + angle = np.arccos(np.dot(vector_to_projection, wall_normal) / ( + np.linalg.norm(vector_to_projection) * np.linalg.norm(wall_normal))) + angle = np.degrees(angle) + + # 如果夹角小于0，将该投影点标记为None + if angle < 90: + intersections[i] = np.array([None, None, None]) + return np.array([]) + + return intersections + + +def is_sunlight_reaching_wall(sun_vec, wall_normal): + # 计算太阳光向量与墙面法线向量之间的夹角 + angle = np.arccos(np.dot(sun_vec, wall_normal) / + (np.linalg.norm(sun_vec) * np.linalg.norm(wall_normal))) + angle = np.degrees(angle) + return angle > 90 + + +def sun_light_vector(azimuth, altitude): + # 将角度转换为弧度 + azimuth_rad = -azimuth+np.pi/2 + altitude_rad = altitude + + # 计算球坐标系中的向量分量 + x = np.cos(altitude_rad) * np.cos(azimuth_rad) + y = np.cos(altitude_rad) * np.sin(azimuth_rad) + z = np.sin(altitude_rad) + + # 太阳光的方向是从太阳指向地球，因此需要反转向量 + return np.array([x, y, -z]) + +def convert_shadows_to_lonlat(all_shadows_coords, center_lon, center_lat): + + a = pd.DataFrame(all_shadows_coords).reset_index() + + b = a[a[0].apply(len) > 0].explode(0).reset_index(drop=True) + b[0] = list(aeqd2lonlat_3d(np.array([list(b[0])]), center_lon, center_lat)[0]) + + c = a[a[0].apply(len) == 0] + + b = b.groupby('index').apply(lambda x: list(x[0])).reset_index() + a = pd.concat([b, c]).sort_values('index').reset_index(drop=True) + return list(a[0]) + +def projections_from_wall_to_wall(merged_data): + + walls_target_shadow = merged_data[merged_data['face'] == False] + + walls_target_shadow = walls_target_shadow[[ + 'building_id_left', 'target_wall_id', 'target_wall', 'date']] + walls_target_shadow['shadow_projections'] = walls_target_shadow['target_wall'] + + # walls_date_shadow = pd.concat([walls_date_shadow, walls_target_shadow]) + + merged_data = merged_data[merged_data['face']] + + def calculate_projection(row): + + return projection_on_wall_single( + row['target_wall_plane'], + row['sun_vector'], + row['shadow_wall'], + row['target_wall_vector'] + ) + + +# 应用修改后的函数计算投影 + merged_data['shadow_projections'] = merged_data.apply( + calculate_projection, axis=1) + merged_data = merged_data.dropna(subset=['shadow_projections']) + +# 或者，如果空的 shadow_projections 以空列表或其他形式存在，您可以这样做： + merged_data = merged_data[merged_data['shadow_projections'].apply( + lambda x: x is not None and len(x) > 0)] + + condition = merged_data.apply(lambda row: + (max([point[0] for point in row['target_wall']]) > min([point[0] for point in row['shadow_projections']])) and + (max([point[2] for point in row['target_wall']]) > min( + [point[2] for point in row['shadow_projections']])), + axis=1) + + merged_data = merged_data[condition] + merged_data = merged_data[[ + 'building_id_left', 'target_wall', 'target_wall_id', 'date', 'shadow_projections']] + + merged_data = pd.concat([merged_data, walls_target_shadow]) + + def aggregate_shadows(group): + walls = group['shadow_projections'].tolist() + aggregated_shadow = cal_multiple_wall_union(walls) + # 返回一个包含所有必要数据的 DataFrame + return pd.DataFrame({ + 'building_id_left': [group.name[0]], + 'target_wall_id': [group.name[1]], + 'date': [group.name[2]], + 'aggregated_shadows': [aggregated_shadow], + 'target_wall': [group['target_wall'].iloc[0]] # 假设仅返回第一个目标墙 + }) + +# 应用 aggregate_shadows 函数 + aggregated_shadows = merged_data.groupby(['building_id_left', 'target_wall_id', 'date']).apply( + aggregate_shadows).reset_index(drop=True) + + +# 展开 DataFrame，因为 apply 返回的是 DataFrame 的列表 + if isinstance(aggregated_shadows.columns, pd.MultiIndex): + aggregated_shadows.columns = [ + '_'.join(col).strip() for col in aggregated_shadows.columns.values] + +# 计算交集并处理结果 + def calculate_intersection(row): + return cal_wall_overlap(row['target_wall'], row['aggregated_shadows'], method='intersection') + + aggregated_shadows['intersection_shadow'] = aggregated_shadows.apply( + calculate_intersection, axis=1) + aggregated_shadows = aggregated_shadows.dropna( + subset=['intersection_shadow']) + +# 删除不需要的列 + result_gdf = aggregated_shadows.drop( + ['aggregated_shadows', 'target_wall'], axis=1) + + result_gdf = result_gdf.dropna(subset=['intersection_shadow']) + +# 或者，如果空的 shadow_projections 以空列表或其他形式存在，您可以这样做： + result_gdf = result_gdf[result_gdf['intersection_shadow'].apply( + lambda x: x is not None and len(x) > 0)] + + return result_gdf + + +def calculate_buildings_shadow_overlap(buildings_gdf, date, precision=3600, padding=1800): + + + buildings_gdf['geometry'] = buildings_gdf['geometry'].apply(make_clockwise) + original_indices = buildings_gdf.index.copy() + + # 确保建筑物数据使用正确的CRS + center_lon, center_lat = buildings_gdf.unary_union.centroid.x, buildings_gdf.unary_union.centroid.y + date_times = get_timetable(center_lon, center_lat, dates=[ + date], precision=precision, padding=padding) + + # 转换建筑物坐标为 AEQD 坐标 + + buildings_aeqd_gdf = buildings_gdf.copy() + buildings_gdf_overlap = buildings_gdf[[ + 'building_id', 'geometry', 'height']] + + for idx, row in buildings_gdf.iterrows(): + if isinstance(row.geometry, Polygon): + lonlat_coords = np.array( + row.geometry.exterior.coords).reshape(1, -1, 2) + elif isinstance(row.geometry, MultiPolygon): + lonlat_coords = np.concatenate([np.array(poly.exterior.coords).reshape( + 1, -1, 2) for poly in row.geometry.geoms], axis=1) + else: + continue + + aeqd_coords = lonlat2aeqd(lonlat_coords, center_lon, center_lat) + buildings_aeqd_gdf.at[idx, 'geometry'] = Polygon( + np.squeeze(aeqd_coords)) + buildings_aeqd_gdf_walls = get_walls(buildings_aeqd_gdf) + buildings_aeqd_gdf_walls = buildings_aeqd_gdf_walls[[ + 'building_id', 'geometry', 'height', 'wall_id']] + walls_shadow = buildings_aeqd_gdf_walls.copy() + + walls_shadow.columns = ['building_id_right', + 'shadow_wall', 'height', 'shadow_wall_id'] + + walls_target = buildings_aeqd_gdf_walls.copy() + + walls_target.columns = ['building_id_left', + 'target_wall', 'height', 'target_wall_id'] + + def polygon_to_points(polygon): + # 提取 Polygon 对象的外部轮廓坐标 + return list(polygon.exterior.coords) + +# 应用转换函数 + walls_target['target_wall'] = walls_target['target_wall'].apply( + polygon_to_points) + + walls_shadow['shadow_wall'] = walls_shadow['shadow_wall'].apply( + polygon_to_points) + walls_target['target_wall_vector'] = walls_target['target_wall'].apply( + calculate_wall_normal_vector) + walls_target['target_wall_plane'] = walls_target['target_wall'].apply( + calculate_wall_plane) + + date_times['date'] = pd.to_datetime(date_times['date']) + merged_data = pd.DataFrame() + + for date_time in date_times['date']: + + sun_position = suncalc.get_position(date_time, center_lon, center_lat) + sun_azimuth = sun_position['azimuth'] + sun_altitude = sun_position['altitude'] + + # 计算所有建筑物的阴影 + shadows_gdf = bdshadow_sunlight(buildings_gdf, date_time) + shadows_gdf.index = original_indices + + # 确保阴影数据也使用相同的CRS + if shadows_gdf.crs is None: + shadows_gdf.set_crs(buildings_gdf.crs, inplace=True) + shadows_gdf = shadows_gdf[['building_id', 'geometry']] + + overlapping = gpd.sjoin(buildings_gdf_overlap, + shadows_gdf, how='left', op='intersects') + + sun_vec = sun_light_vector(sun_azimuth, sun_altitude) + walls_target['sun_vector'] = [sun_vec] * len(walls_target) + + def have_sun(row): + + return is_sunlight_reaching_wall( + row['sun_vector'], + row['target_wall_vector'] + ) + + walls_target['face'] = walls_target.apply(have_sun, axis=1) + + walls_target['date'] = date_time + overlapping = pd.merge( + overlapping, walls_target, on='building_id_left') + + overlapping = pd.merge(overlapping, walls_shadow, + on='building_id_right') + + overlapping = overlapping[-((overlapping['building_id_left'] == overlapping['building_id_right']) & + (overlapping['target_wall_id'] == overlapping['shadow_wall_id']))] + + merged_data = pd.concat([merged_data, overlapping]) + # walls_target = walls_target.drop(['date','face','sun_vector'],axis=1) + + final_merged_data = projections_from_wall_to_wall(merged_data) + + final_merged_data['intersection_shadow_lonlat'] = convert_shadows_to_lonlat( + final_merged_data['intersection_shadow'].tolist(), center_lon, center_lat) + final_merged_data = final_merged_data.dropna( + subset=['intersection_shadow_lonlat']) + + def create_polygon_from_coords(coords): + # 仅当坐标列表不为空时创建 Polygon + # if coords is not None and len(coords) > 2: + return Polygon(coords) + # return None + +# 将 intersection_shadow 转换为 Polygon 对象 + final_merged_data['intersection_shadow_polygon'] = final_merged_data['intersection_shadow_lonlat'].apply( + create_polygon_from_coords) + # final_merged_data['target_wall'] = final_merged_data['target_wall'].apply(create_polygon_from_coords) + final_merged_data = final_merged_data.rename( + columns={'building_id_left': 'building_index'}) + # final_merged_data = final_merged_data.dropna(subset=['intersection_shadow_polygon']) + # all_shadows_gdf = pd.concat([all_shadows_gdf, final_merged_data]) + + final_merged_data = final_merged_data.drop( + ['intersection_shadow', 'intersection_shadow_lonlat'], axis=1) + + return final_merged_data + +def cal_sunshine_facade(buildings_gdf, day, precision=3600, padding=1800): + + # 计算阴影重叠情况 + final_shadow = calculate_buildings_shadow_overlap( + buildings_gdf, day, precision=precision, padding=padding) + + final_shadow['building_id'] = final_shadow['building_index'] + + # 求最大光照时长 + lon, lat = buildings_gdf['geometry'].iloc[0].bounds[:2] + date = pd.to_datetime(day+' 12:45:33.959797119') + times = get_times(date, lon, lat) + date_sunrise = times['sunrise'] + data_sunset = times['sunset'] + timestamp_sunrise = pd.Series(date_sunrise).astype('int') + timestamp_sunset = pd.Series(data_sunset).astype('int') + sunlighthour = ( + timestamp_sunset.iloc[0]-timestamp_sunrise.iloc[0])/(1000000000*3600) + + # 从阴影重叠情况计算光照时长 + final_shadows_oneday = final_shadow.copy() + buildings_walls = get_walls(buildings_gdf) + buildings_walls = buildings_walls.rename( + columns={'wall_id': 'target_wall_id'}) + buildings_walls = buildings_walls[[ + 'building_id', 'target_wall_id', 'geometry']] + final_shadows_oneday = final_shadows_oneday.rename( + columns={'intersection_shadow_polygon': 'geometry'}) + final_shadows_oneday = pd.concat([final_shadows_oneday, buildings_walls]) + final_shadows_oneday = final_shadows_oneday[final_shadows_oneday['geometry'].apply( + lambda x: has_normal(x.exterior.coords))] + final_shadows_sunshinetime = final_shadows_oneday.groupby( + ['building_id', 'target_wall_id']).apply(cal_multiple_wall_overlap_count) + + final_shadows_sunshinetime = final_shadows_sunshinetime.reset_index() + + final_shadows_sunshinetime['time'] = ( + final_shadows_sunshinetime['count']-1)*precision + final_shadows_sunshinetime['Hour'] = sunlighthour - \ + final_shadows_sunshinetime['time']/3600 + final_shadows_sunshinetime.loc[final_shadows_sunshinetime['Hour'] + <= 0, 'Hour'] = 0 + return final_shadows_sunshinetime + + +def get_walls(buildings_gdf): + # 把建筑物的墙面拆分成单独的面 + def get_wall(r): + wall_coords = list(r['geometry'].exterior.coords) + walls = [] + for i in range(len(wall_coords)-1): + walls.append([Polygon([list(wall_coords[i])+[0], + list(wall_coords[i+1])+[0], + list(wall_coords[i+1])+[r['height']], + list(wall_coords[i])+[r['height']], + ]), i]) + return walls + buildings_gdf['walls'] = buildings_gdf.apply(get_wall, axis=1) + + buildings_walls = buildings_gdf.explode('walls') + buildings_walls['wall_id'] = buildings_walls['walls'].apply(lambda x: x[1]) + buildings_walls['geometry'] = buildings_walls['walls'].apply( + lambda x: x[0]) + return buildings_walls.drop('walls', axis=1) diff --git a/src/pybdshadow/get_buildings.py b/src/pybdshadow/get_buildings.py new file mode 100644 index 0000000..d93a020 --- /dev/null +++ b/src/pybdshadow/get_buildings.py @@ -0,0 +1,279 @@ +import requests +from vt2geojson.tools import vt_bytes_to_geojson +import pandas as pd +import geopandas as gpd +import transbigdata as tbd +from .preprocess import bd_preprocess +from tqdm import tqdm +import math +from retrying import retry +from requests.exceptions import RequestException + +def deg2num(lat_deg, lon_deg, zoom): + ''' + Calculate xy tiles from coordinates + + Parameters + ------- + lon_deg : number + Longitude + lat_deg : number + Latitude + zoom : Int + Zoom level of the map + ''' + lat_rad = math.radians(lat_deg) + n = 2.0 ** zoom + xtile = int((lon_deg + 180.0) / 360.0 * n) + ytile = int((1.0 - math.log(math.tan(lat_rad) + + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n) + return (xtile, ytile) + + + +def is_request_exception(e): + return issubclass(type(e),RequestException) + +@retry(retry_on_exception=is_request_exception,wrap_exception=False, stop_max_attempt_number=300) +def safe_request(url, **kwargs): + return requests.get(url, **kwargs) + +def getbd(x,y,z,MAPBOX_ACCESS_TOKEN): + ''' + Get buildings from mapbox vector tiles + + Parameters + ------- + x : Int + x tile number + y : Int + y tile number + z : Int + zoom level of the map + MAPBOX_ACCESS_TOKEN : str + Mapbox access token + + Return + ---------- + building : GeoDataFrame + buildings in the tile + ''' + try: + url = f"https://api.mapbox.com/v4/mapbox.mapbox-streets-v8,mapbox.mapbox-terrain-v2,mapbox.mapbox-bathymetry-v2/{z}/{x}/{y}.vector.pbf?sku=101vMyxQx9v3Q&access_token={MAPBOX_ACCESS_TOKEN}" + + r = safe_request(url, timeout=10) + assert r.status_code == 200, r.content + vt_content = r.content + features = vt_bytes_to_geojson(vt_content, x, y, z) + gdf = gpd.GeoDataFrame.from_features(features) + building = gdf[gdf['height']>0][['geometry', 'height','type']] + except: + building = pd.DataFrame() + return building + + +def get_tiles_by_lonlat(lon1,lat1,lon2,lat2,z): + ''' + Get tiles by lonlat + + Parameters + ------- + lon1 : number + Longitude of the first point + lat1 : number + Latitude of the first point + lon2 : number + Longitude of the second point + lat2 : number + Latitude of the second point + z : Int + Zoom level of the map + + Return + ---------- + tiles : DataFrame + Tiles in the area + ''' + x1,y1 = deg2num(lat1, lon1, z) + x2,y2 = deg2num(lat2, lon2, z) + x_min = min(x1,x2) + x_max = max(x1,x2) + y_min = min(y1,y2) + y_max = max(y1,y2) + tiles = pd.DataFrame(range(x_min,x_max+1), columns=['x']).assign(foo=1).merge(pd.DataFrame(range(y_min,y_max+1), columns=['y']).assign(foo=1)).drop('foo', axis=1).assign(z=z) + return tiles + +def get_tiles_by_polygon(polygon,z): + ''' + Get tiles by polygon + + Parameters + ------- + polygon : GeoDataFrame of Polygon or MultiPolygon + Polygon of the area + z : Int + Zoom level of the map + + Return + ---------- + tiles : DataFrame + Tiles in the area + ''' + grid,params = tbd.area_to_grid(polygon,accuracy=400) + grid['lon'] = grid.centroid.x + grid['lat'] = grid.centroid.y + a = grid.apply(lambda x: deg2num(x.lat, x.lon, z), axis=1) + grid['x'] = a.apply(lambda a:a[0]) + grid['y'] = a.apply(lambda a:a[1]) + grid['z'] = z + tiles = grid[['x','y','z']].drop_duplicates() + return tiles + +def get_buildings_threading(tiles,MAPBOX_ACCESS_TOKEN,merge=False,num_threads=100): + ''' + Get buildings by threading + + Parameters + ------- + tiles : DataFrame + Tiles in the area + MAPBOX_ACCESS_TOKEN : str + Mapbox access token + merge : bool + whether to merge buildings in the same grid + num_threads : Int + number of threads + + Return + ---------- + building : GeoDataFrame + buildings in the area + ''' + def merge_building(building): + building = building.groupby(['height','type']).apply(lambda r:r.unary_union).reset_index() + building.columns = ['height','type','geometry'] + building = gpd.GeoDataFrame(building,geometry = 'geometry') + building = bd_preprocess(building) + return building + + # 这是修改后的 getbd_tojson 函数 + def getbd_tojson(data, MAPBOX_ACCESS_TOKEN, pbar, results): + for j in range(len(data)): + r = data.iloc[j] + x, y, z = r['x'], r['y'], r['z'] + try: + url = f"https://api.mapbox.com/v4/mapbox.mapbox-streets-v8,mapbox.mapbox-terrain-v2,mapbox.mapbox-bathymetry-v2/{z}/{x}/{y}.vector.pbf?sku=101vMyxQx9v3Q&access_token={MAPBOX_ACCESS_TOKEN}" + r = safe_request(url, timeout=10) + assert r.status_code == 200, r.content + vt_content = r.content + features = vt_bytes_to_geojson(vt_content, x, y, z) + gdf = gpd.GeoDataFrame.from_features(features) + building = gdf[gdf['height'] > 0][['geometry', 'height', 'type']] + results.append(building) # 将结果添加到全局列表 + except: + pass + finally: + pbar.update() + + # 主程序 + import threading + import os + # 主程序 + # 分割数据 + grid = tiles.copy() + bins = num_threads + + grid['tmpid'] = range(len(grid)) + grid['group_num'] = pd.cut(grid['tmpid'], bins, precision=2, labels=range(bins)) + + # 创建进度条 + pbar = tqdm(total=len(grid), desc='Downloading Buildings: ') + + # 存储结果的全局列表 + results = [] + + # 划分线程 + threads = [] + for i in range(bins): + data = grid[grid['group_num'] == i] + threads.append(threading.Thread(target=getbd_tojson, args=(data, MAPBOX_ACCESS_TOKEN, pbar, results))) + + # 线程开始 + for t in threads: + t.setDaemon(True) + t.start() + for t in threads: + t.join() + + # 关闭进度条 + pbar.close() + threads.clear() + + # 合并数据 + building = pd.concat(results) + + if merge: + #再做一次聚合，分栅格聚合建筑 + building['x'] = building.centroid.x + building['y'] = building.centroid.y + params = tbd.area_to_params(building['geometry'].iloc[0].bounds) + building['LONCOL'],building['LATCOL'] = tbd.GPS_to_grid(building['x'],building['y'],params) + building['tile'] = building['LONCOL'].astype(str)+'_'+building['LATCOL'].astype(str) + building = building.groupby(['tile','type']).apply(merge_building).reset_index(drop=True) + + building = building[['geometry','height','type']] + building['building_id'] = range(len(building)) + + return building + +def get_buildings_by_bounds(lon1,lat1,lon2,lat2,MAPBOX_ACCESS_TOKEN,merge=False): + ''' + Get buildings by bounds + + Parameters + ------- + lon1 : number + Longitude of the first point + lat1 : number + Latitude of the first point + lon2 : number + Longitude of the second point + lat2 : number + Latitude of the second point + MAPBOX_ACCESS_TOKEN : str + Mapbox access token + merge : bool + whether to merge buildings in the same grid + + Return + ---------- + building : GeoDataFrame + buildings in the area + ''' + tiles = get_tiles_by_lonlat(lon1,lat1,lon2,lat2,16) + building = get_buildings_threading(tiles,MAPBOX_ACCESS_TOKEN,merge) + building = bd_preprocess(building) + return building + +def get_buildings_by_polygon(polygon,MAPBOX_ACCESS_TOKEN,merge=False): + ''' + Get buildings by polygon + + Parameters + ------- + polygon : GeoDataFrame of Polygon or MultiPolygon + Polygon of the area + MAPBOX_ACCESS_TOKEN : str + Mapbox access token + merge : bool + whether to merge buildings in the same grid + + Return + ---------- + building : GeoDataFrame + buildings in the area + ''' + tiles = get_tiles_by_polygon(polygon,16) + building = get_buildings_threading(tiles,MAPBOX_ACCESS_TOKEN,merge) + building = bd_preprocess(building) + return building \ No newline at end of file diff --git a/src/pybdshadow/preprocess.py b/src/pybdshadow/preprocess.py index 487f8af..9ad5165 100644 --- a/src/pybdshadow/preprocess.py +++ b/src/pybdshadow/preprocess.py @@ -32,7 +32,7 @@ import shapely import pandas as pd import geopandas as gpd - +from shapely.geometry import MultiPolygon def bd_preprocess(buildings, height=''): ''' @@ -79,6 +79,7 @@ def bd_preprocess(buildings, height=''): allbds['geometry'] = allbds.buffer(0) else: allbds = gpd.GeoDataFrame() + allbds.crs = 'EPSG:4326' return allbds def gdf_difference(gdf_a,gdf_b,col = 'building_id'): @@ -89,7 +90,7 @@ def gdf_difference(gdf_a,gdf_b,col = 'building_id'): gdfb = gdf_b.copy() gdfb = gdfb[['geometry']] #判断重叠 - from shapely.geometry import MultiPolygon + gdfa.crs = gdfb.crs gdfb = gpd.sjoin(gdfb,gdfa).groupby([col])['geometry'].apply( lambda df: MultiPolygon(list(df)).buffer(0)).reset_index() @@ -116,7 +117,6 @@ def gdf_intersect(gdf_a,gdf_b,col = 'building_id'): gdfb = gdf_b.copy() gdfb = gdfb[['geometry']] #判断重叠 - from shapely.geometry import MultiPolygon gdfa.crs = gdfb.crs gdfb = gpd.sjoin(gdfb,gdfa).groupby([col])['geometry'].apply( lambda df: MultiPolygon(list(df)).buffer(0)).reset_index() diff --git a/src/pybdshadow/pybdshadow.py b/src/pybdshadow/pybdshadow.py index f53e1a6..9a8588f 100644 --- a/src/pybdshadow/pybdshadow.py +++ b/src/pybdshadow/pybdshadow.py @@ -32,7 +32,7 @@ import pandas as pd import geopandas as gpd from suncalc import get_position -from shapely.geometry import Polygon,LineString, MultiPolygon +from shapely.geometry import Polygon, MultiPolygon import math import numpy as np from .utils import ( @@ -63,7 +63,7 @@ def calSunShadow_vector(shape, shapeHeight, sunPosition): # transform coordinate system meanlon = shape[:,:,0].mean() meanlat = shape[:,:,1].mean() - shape = lonlat2aeqd(shape) + shape = lonlat2aeqd(shape,meanlon,meanlat) azimuth = sunPosition['azimuth'] altitude = sunPosition['altitude'] diff --git a/src/pybdshadow/utils.py b/src/pybdshadow/utils.py index 1d634d2..f993b5c 100644 --- a/src/pybdshadow/utils.py +++ b/src/pybdshadow/utils.py @@ -31,8 +31,24 @@ """ import numpy as np - -def lonlat2aeqd(lonlat): +import shapely +import geopandas as gpd +import pandas as pd +from pyproj import CRS,Transformer +from shapely.geometry import Polygon + +def extrude_poly(poly,h): + poly_coords = np.array(poly.exterior.coords) + poly_coords = np.c_[poly_coords, np.ones(poly_coords.shape[0])*h] + return Polygon(poly_coords) + +def make_clockwise(polygon): + if polygon.exterior.is_ccw: + return polygon + else: + return Polygon(list(polygon.exterior.coords)[::-1]) + +def lonlat2aeqd(lonlat, center_lon, center_lat): ''' Convert longitude and latitude to azimuthal equidistant projection coordinates. @@ -58,16 +74,31 @@ def lonlat2aeqd(lonlat): [[-48243.5939812 , -55322.02388971], [ 47752.57582735, 55538.86412435]]]) ''' - meanlon = lonlat[:,:,0].mean() - meanlat = lonlat[:,:,1].mean() - from pyproj import CRS - epsg = CRS.from_proj4("+proj=aeqd +lat_0="+str(meanlat)+" +lon_0="+str(meanlon)+" +datum=WGS84") - from pyproj import Transformer - transformer = Transformer.from_crs("EPSG:4326", epsg,always_xy = True) - proj_coords = transformer.transform(lonlat[:,:,0], lonlat[:,:,1]) - proj_coords = np.array(proj_coords).transpose([1,2,0]) + epsg = CRS.from_proj4("+proj=aeqd +lat_0="+str(center_lat) + + " +lon_0="+str(center_lon)+" +datum=WGS84") + transformer = Transformer.from_crs("EPSG:4326", epsg, always_xy=True) + proj_coords = transformer.transform(lonlat[:, :, 0], lonlat[:, :, 1]) + proj_coords = np.array(proj_coords).transpose([1, 2, 0]) return proj_coords +def aeqd2lonlat_3d(proj_coords, meanlon, meanlat): + + # 提取 xy 坐标和 z 坐标 + xy_coords = proj_coords[:, :, :2] + z_coords = proj_coords[:, :, 2] if proj_coords.shape[2] > 2 else np.zeros( + xy_coords.shape[:2]) + + # 定义转换器 + epsg = CRS.from_proj4("+proj=aeqd +lat_0=" + str(meanlat) + + " +lon_0=" + str(meanlon) + " +datum=WGS84") + transformer = Transformer.from_crs(epsg, "EPSG:4326", always_xy=True) + + # 转换 xy 坐标 + lon, lat = transformer.transform(xy_coords[:, :, 0], xy_coords[:, :, 1]) + + # 将转换后的坐标和原始 z 坐标组合 + lonlat = np.dstack([lon, lat, z_coords]) + return lonlat def aeqd2lonlat(proj_coords,meanlon,meanlat): ''' @@ -103,11 +134,68 @@ def aeqd2lonlat(proj_coords,meanlon,meanlat): [[120., 30.], [121., 31.]]]) ''' - from pyproj import CRS + epsg = CRS.from_proj4("+proj=aeqd +lat_0="+str(meanlat)+" +lon_0="+str(meanlon)+" +datum=WGS84") - from pyproj import Transformer transformer = Transformer.from_crs( epsg,"EPSG:4326",always_xy = True) lonlat = transformer.transform(proj_coords[:,:,0], proj_coords[:,:,1]) lonlat = np.array(lonlat).transpose([1,2,0]) return lonlat +def calculate_normal(points): + points = np.array(points) + if points.shape[0] < 3: + raise ValueError("墙至少需要三个点。") + + for i in range(points.shape[0]): + for j in range(i + 1, points.shape[0]): + for k in range(j + 1, points.shape[0]): + vector1 = points[j] - points[i] + vector2 = points[k] - points[i] + normal = np.cross(vector1, vector2) + if np.linalg.norm(normal) != 0: + return normal / np.linalg.norm(normal) + + raise ValueError("该墙所有点共线，无法计算法向量。") + +def has_normal(points): + # 将点列表转换为NumPy数组以便处理 + points = np.array(points) + + # 需要至少三个点来形成一个平面 + if points.shape[0] < 3: + return False + + # 寻找不共线的三个点 + for i in range(points.shape[0]): + for j in range(i+1, points.shape[0]): + for k in range(j+1, points.shape[0]): + # 计算两个向量 + vector1 = points[j] - points[i] + vector2 = points[k] - points[i] + + # 计算叉乘 + normal = np.cross(vector1, vector2) + + # 检查法向量是否非零（即点不共线） + if np.linalg.norm(normal) != 0: + # 返回归一化的法向量 + return True + return False + +def count_overlapping_features(gdf,buffer = True): + # 计算多边形的重叠次数 + if buffer: + bounds = gdf.geometry.buffer(1e-9).exterior.unary_union + else: + bounds = gdf.geometry.exterior.unary_union + + new_polys = list(shapely.ops.polygonize(bounds)) + new_gdf = gpd.GeoDataFrame(geometry=new_polys) + new_gdf['id'] = range(len(new_gdf)) + new_gdf_centroid = new_gdf.copy() + new_gdf_centroid['geometry'] = new_gdf.geometry.representative_point() + overlapcount = gpd.sjoin(new_gdf_centroid, gdf) + overlapcount = overlapcount.groupby( + ['id'])['index_right'].count().rename('count').reset_index() + out_gdf = pd.merge(new_gdf, overlapcount) + return out_gdf \ No newline at end of file diff --git a/src/pybdshadow/visualization.py b/src/pybdshadow/visualization.py index 22ef4b7..5a1336e 100644 --- a/src/pybdshadow/visualization.py +++ b/src/pybdshadow/visualization.py @@ -32,13 +32,15 @@ import numpy as np import geopandas as gpd +from shapely.geometry import Polygon def show_bdshadow(buildings=gpd.GeoDataFrame(), shadows=gpd.GeoDataFrame(), ad=gpd.GeoDataFrame(), ad_visualArea=gpd.GeoDataFrame(), height='height', - zoom='auto'): + zoom='auto', + vis_height = 800): ''' Visualize the building and shadow with keplergl. @@ -391,5 +393,144 @@ def show_bdshadow(buildings=gpd.GeoDataFrame(), 'building': True, 'water': True, 'land': True}, - 'mapStyles': {}}}}, data=vmapdata, height=500) + 'mapStyles': {}}}}, data=vmapdata, height=vis_height) return vmap + + + +def show_sunshine(sunshine=gpd.GeoDataFrame(), + zoom='auto',vis_height = 800): + ''' + Visualize the sunshine with keplergl. + + Parameters + -------------------- + sunshine : GeoDataFrame + sunshine. coordinate system should be WGS84 + zoom : number + Zoom level of the map + + Return + -------------------- + vmap : keplergl.keplergl.KeplerGl + Visualizations provided by keplergl + ''' + def offset_wall(wall_poly): + wall_coords = np.array(wall_poly.exterior.coords) + wall_coords[:,0]+=wall_coords[:,2]*0.000000001 + wall_coords[:,1]+=wall_coords[:,2]*0.000000001 + return Polygon(wall_coords) + sunshine = sunshine.copy() + sunshine['geometry'] = sunshine['geometry'].apply(offset_wall) + vmapdata = {} + layers = [] + + bdcentroid = sunshine['geometry'].bounds[[ + 'minx', 'miny', 'maxx', 'maxy']] + lon_center, lat_center = bdcentroid['minx'].mean( + ), bdcentroid['miny'].mean() + lon_min, lon_max = bdcentroid['minx'].min(), bdcentroid['maxx'].max() + vmapdata['sunshine'] = sunshine + + + layers.append( + {'id': 'lz48o4', + 'type': 'geojson', + 'config': { + 'dataId': 'sunshine', + 'label': 'sunshine', + 'color': [73, 73, 73], + 'highlightColor': [252, 242, 26, 255], + 'columns': {'geojson': 'geometry'}, + 'isVisible': True, + 'visConfig': { + 'opacity': 1, + 'strokeOpacity': 1, + 'thickness': 0.5, + 'strokeColor': [255, 153, 31], + 'colorRange': {'name': 'UberPool 9', + 'type': 'sequential', + 'category': 'Uber', + 'colors': ['#2C51BE', + '#482BBD', + '#7A0DA6', + '#AE0E7F', + '#CF1750', + '#E31A1A', + '#FD7900', + '#FAC200', + '#FAE300'], + 'reversed': False}, + 'strokeColorRange': {'name': 'Global Warming', + 'type': 'sequential', + 'category': 'Uber', + 'colors': ['#5A1846', + '#900C3F', + '#C70039', + '#E3611C', + '#F1920E', + '#FFC300']}, + 'radius': 10, + 'sizeRange': [0, 10], + 'radiusRange': [0, 50], + 'heightRange': [0, 500], + 'elevationScale': 5, + 'enableElevationZoomFactor': True, + 'stroked': False, + 'filled': True, + 'enable3d': False, + 'wireframe': False}, + 'hidden': False, + 'textLabel': [{ + 'field': None, + 'color': [255, 255, 255], + 'size': 18, + 'offset': [0, 0], + 'anchor': 'start', + 'alignment': 'center'}]}, + 'visualChannels': { + 'colorField': {'name': 'Hour', 'type': 'real'}, + 'colorScale': 'quantize', + 'strokeColorField': None, + 'strokeColorScale': 'quantize', + 'sizeField': None, + 'sizeScale': 'linear', + 'heightField': None, + 'heightScale': 'linear', + 'radiusField': None, + 'radiusScale': 'linear'}}) + try: + from keplergl import KeplerGl + except ImportError: + raise ImportError( + "Please install keplergl, run " + "the following code in cmd: pip install keplergl") + + if zoom == 'auto': + zoom = 10.5-np.log(lon_max-lon_min)/np.log(2) + vmap = KeplerGl(config={ + 'version': 'v1', + 'config': { + 'visState': { + 'filters': [], + 'layers': layers, + 'layerBlending': 'normal', + 'animationConfig': {'currentTime': None, 'speed': 1}}, + 'mapState': {'bearing': 30, + 'dragRotate': True, + 'latitude': lat_center, + 'longitude': lon_center, + 'pitch': 50, + 'zoom': zoom, + 'isSplit': False}, + 'mapStyle': {'styleType': 'light', + 'topLayerGroups': {}, + 'visibleLayerGroups': {'label': True, + 'road': True, + 'border': False, + 'building': True, + 'water': True, + 'land': True}, + 'mapStyles': {}}}}, data=vmapdata, height=vis_height) + return vmap +