Merge pull request #5 from ni1o1/0.2.3

0.2.3

docs/source/bdshadow.rst

@@ -9,7 +9,7 @@ Building shadow
 Shadow from sunlight
 --------------------------------------
 
-.. function:: pybdshadow.bdshadow_sunlight(buildings, date, merge=False, height='height', ground=0)
+.. function:: pybdshadow.bdshadow_sunlight(buildings, date,  height='height', roof=False,include_building = True,ground=0)
 
 Calculate the sunlight shadow of the buildings.
 
@@ -19,10 +19,12 @@ buildings : GeoDataFrame
     Buildings. coordinate system should be WGS84
 date : datetime
     Datetime
-merge : bool
-    whether to merge the wall shadows into the building shadows
 height : string
     Column name of building height
+roof : bool
+    whether to calculate the roof shadows
+include_building : bool
+    whether the shadow include building outline
 ground : number
     Height of the ground
 

docs/source/conf.py

@@ -22,8 +22,8 @@
 author = 'Qing Yu'
 
 # The full version, including alpha/beta/rc tags
-release = '0.2.2'
-version = '0.2.2'
+release = '0.2.3'
+version = '0.2.3'
 html_logo = "_static/logo-wordmark-light.png"
 html_favicon = '_static/logo.ico'
 # -- General configuration ---------------------------------------------------

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="pybdshadow",
-    version="0.2.2",
+    version="0.2.3",
     author="Qing Yu",
     author_email="[email protected]",
     description="Python package to generate building shadow geometry",

src/pybdshadow/__init__.py

@@ -32,7 +32,7 @@
 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 """
 
-__version__ = '0.2.2'
+__version__ = '0.2.3'
 __author__ = 'Qing Yu <[email protected]>'
 
 # module level doc-string

src/pybdshadow/preprocess.py

@@ -71,3 +71,55 @@ def bd_preprocess(buildings, height='height'):
     allbds = pd.concat(allbds)
     allbds['building_id'] = range(len(allbds))
     return allbds
+
+def gdf_difference(gdf_a,gdf_b,col = 'building_id'):
+    '''
+    difference gdf_b from gdf_a
+    '''
+    gdfa = gdf_a.copy()
+    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()
+    #分割有重叠和无重叠的
+    gdfb['tmp'] = 1
+    gdfa_1 = pd.merge(gdfa,gdfb[[col,'tmp']],how = 'left')
+    gdfa = gdfa_1[gdfa_1['tmp'] == 1].drop('tmp',axis = 1)
+    gdfa_notintersected = gdfa_1[gdfa_1['tmp'].isnull()].drop('tmp',axis = 1)
+    #对有重叠的进行裁剪
+    gdfa = gdfa.sort_values(by = col).set_index(col)
+    gdfb = gdfb.sort_values(by = col).set_index(col)
+    gdfa.crs = gdfb.crs
+    gdfa['geometry'] = gdfa.difference(gdfb)
+    gdfa = gdfa.reset_index()
+    #拼合
+    gdfa = pd.concat([gdfa,gdfa_notintersected])
+    return gdfa
+
+def gdf_intersect(gdf_a,gdf_b,col = 'building_id'):
+    '''
+    intersect gdf_b from gdf_a
+    '''
+    gdfa = gdf_a.copy()
+    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()
+    #分割有重叠和无重叠的
+    gdfb['tmp'] = 1
+    gdfa_1 = pd.merge(gdfa,gdfb[[col,'tmp']],how = 'left')
+    gdfa = gdfa_1[gdfa_1['tmp'] == 1].drop('tmp',axis = 1)
+    #对有重叠的进行裁剪
+    gdfa = gdfa.sort_values(by = col).set_index(col)
+    gdfb = gdfb.sort_values(by = col).set_index(col)
+    gdfa.crs = gdfb.crs
+    gdfa['geometry'] = gdfa.intersection(gdfb)
+    gdfa = gdfa.reset_index()
+
+    return gdfa

src/pybdshadow/pybdshadow.py

@@ -39,6 +39,7 @@
     lonlat_mercator_vector,
     mercator_lonlat_vector
 )
+from .preprocess import gdf_difference,gdf_intersect
 
 
 def calSunShadow_vector(shape, shapeHeight, sunPosition):
@@ -76,23 +77,27 @@ def calSunShadow_vector(shape, shapeHeight, sunPosition):
     return shadowShape
 
 
-def bdshadow_sunlight(buildings, date, merge=True, height='height', ground=0):
+def bdshadow_sunlight(buildings, date,  height='height', roof=False,include_building = True,ground=0):
     '''
     Calculate the sunlight shadow of the buildings.
 
     **Parameters**
+
     buildings : GeoDataFrame
         Buildings. coordinate system should be WGS84
     date : datetime
         Datetime
-    merge : bool
-        whether to merge the wall shadows into the building shadows
     height : string
         Column name of building height
+    roof : bool
+        whether to calculate the roof shadows
+    include_building : bool
+        whether the shadow include building outline
     ground : number
         Height of the ground
 
     **Return**
+
     shadows : GeoDataFrame
         Building shadow
     '''
@@ -126,21 +131,99 @@ def bdshadow_sunlight(buildings, date, merge=True, height='height', ground=0):
     walls = walls[['x1', 'y1', 'x2', 'y2', 'building_id', 'height']]
     walls['wall'] = walls.apply(lambda r: [[r['x1'], r['y1']],
                                            [r['x2'], r['y2']]], axis=1)
-    walls_shape = np.array(list(walls['wall']))
+
+    ground_shadow = walls.copy()
+    walls_shape = np.array(list(ground_shadow['wall']))
 
     # calculate shadow for walls
     shadowShape = calSunShadow_vector(
-        walls_shape, walls['height'].values, sunPosition)
-
-    walls['geometry'] = list(shadowShape)
-    walls['geometry'] = walls['geometry'].apply(lambda r: Polygon(r))
-    walls = gpd.GeoDataFrame(walls)
-    walls = pd.concat([walls, building])
-    if merge:
-        walls = walls.groupby(['building_id'])['geometry'].apply(
-            lambda df: MultiPolygon(list(df)).buffer(0)).reset_index()
-
-    return walls
+        walls_shape, ground_shadow['height'].values, sunPosition)
+
+    ground_shadow['geometry'] = list(shadowShape)
+    ground_shadow['geometry'] = ground_shadow['geometry'].apply(
+        lambda r: Polygon(r))
+    ground_shadow = gpd.GeoDataFrame(ground_shadow)
+
+
+
+    ground_shadow = pd.concat([ground_shadow, building])
+    ground_shadow = ground_shadow.groupby(['building_id'])['geometry'].apply(
+        lambda df: MultiPolygon(list(df)).buffer(0)).reset_index()
+
+    ground_shadow['height'] = 0
+    ground_shadow['type'] = 'ground'
+
+    if not roof:
+        if not include_building:
+            #从地面阴影裁剪建筑轮廓
+            ground_shadow = gdf_difference(ground_shadow,buildings)
+        return ground_shadow
+    else:
+        def calwall_shadow(walldata, building):
+            walls = walldata.copy()
+            walls_shape = np.array(list(walls['wall']))
+            # calculate shadow for walls
+            shadowShape = calSunShadow_vector(
+                walls_shape, walls['height'].values, sunPosition)
+            walls['geometry'] = list(shadowShape)
+            walls['geometry'] = walls['geometry'].apply(lambda r: Polygon(r))
+            walls = gpd.GeoDataFrame(walls)
+            walls = pd.concat([walls, building])
+
+            walls = walls.groupby(['building_id'])['geometry'].apply(
+                lambda df: MultiPolygon(list(df)).buffer(0)).reset_index()
+            return walls
+
+        # 计算屋顶阴影
+        roof_shadows = []
+        for roof_height in walls[height].drop_duplicates():
+            # 高于给定高度的墙
+            walls_high = walls[walls[height] > roof_height].copy()
+            if len(walls_high) == 0:
+                continue
+            walls_high[height] -= roof_height
+            # 高于给定高度的建筑
+            building_high = building[building[height] > roof_height].copy()
+            if len(building_high) == 0:
+                continue
+            building_high[height] -= roof_height
+            # 所有建筑在此高度的阴影
+            building_shadow_height = calwall_shadow(walls_high, building_high)
+            # 在此高度的建筑屋顶
+            building_roof = building[building[height] == roof_height].copy()
+            building_shadow_height.crs = building_roof.crs
+            # 取有遮挡的阴影
+            building_shadow_height = gpd.sjoin(
+                building_shadow_height, building_roof)
+            if len(building_shadow_height) == 0:
+                continue
+            # 与屋顶做交集
+            building_roof = gdf_intersect(building_roof,building_shadow_height)
+
+            # 再减去这个高度以上的建筑
+            building_higher = building[building[height] > roof_height].copy()
+            building_roof = gdf_difference(building_roof,building_higher)
+
+            #给出高度信息
+            building_roof['height'] = roof_height
+            building_roof = building_roof[-building_roof['geometry'].is_empty]
+
+            roof_shadows.append(building_roof)
+        if len(roof_shadows) == 0:
+            roof_shadow = gpd.GeoDataFrame()
+        else:
+            roof_shadow = pd.concat(roof_shadows)[
+                ['height', 'building_id', 'geometry']]
+            roof_shadow['type'] = 'roof'
+
+        if not include_building:
+            #从地面阴影裁剪建筑轮廓
+            ground_shadow = gdf_difference(ground_shadow,buildings)
+
+        shadows = pd.concat([roof_shadow, ground_shadow])
+        shadows.crs = None
+        shadows['geometry'] = shadows.buffer(0.000001).buffer(-0.000001)
+        return shadows
 
 
 def calPointLightShadow_vector(shape, shapeHeight, pointLight):
@@ -183,6 +266,7 @@ def bdshadow_pointlight(buildings,
     Calculate the sunlight shadow of the buildings.
 
     **Parameters**
+    
     buildings : GeoDataFrame
         Buildings. coordinate system should be WGS84
     pointlon,pointlat,pointheight : float
@@ -197,6 +281,7 @@ def bdshadow_pointlight(buildings,
         Height of the ground
 
     **Return**
+
     shadows : GeoDataFrame
         Building shadow
     '''

src/pybdshadow/tests/test_pybdshadow.py

@@ -7,43 +7,49 @@
 
 class Testpybdshadow:
     def test_bdshadow_sunlight(self):
+
         buildings = gpd.GeoDataFrame({
-            'height': [42],
+            'height': [42, 9],
             'geometry': [
                 Polygon([(139.698311, 35.533796),
-                         (139.698311,
-                          35.533642),
-                         (139.699075,
-                          35.533637),
-                         (139.699079,
-                          35.53417),
-                         (139.698891,
-                          35.53417),
-                         (139.698888,
-                          35.533794),
-                         (139.698311, 35.533796)])]})
-        date = pd.to_datetime('2015-01-01 02:45:33.959797119')
+                        (139.698311,
+                            35.533642),
+                        (139.699075,
+                            35.533637),
+                        (139.699079,
+                            35.53417),
+                        (139.698891,
+                            35.53417),
+                        (139.698888,
+                            35.533794),
+                        (139.698311, 35.533796)]),
+                Polygon([(139.69799, 35.534175),
+                        (139.697988, 35.53389),
+                        (139.698814, 35.533885),
+                        (139.698816, 35.534171),
+                        (139.69799, 35.534175)])]})
+        date = pd.to_datetime('2015-01-01 03:45:33.959797119')
 
         buildings = pybdshadow.bd_preprocess(buildings)
 
         buildingshadow = pybdshadow.bdshadow_sunlight(
-            buildings, date, merge=True)
+            buildings, date, roof=True, include_building=False)
 
         area = buildingshadow['geometry'].iloc[0]
         area = np.array(area.exterior.coords)
-        truth = np.array([[139.698311,  35.533796],
-                          [139.69831102,  35.533796],
-                          [139.69831102,  35.533796],
-                          [139.69831239,  35.53429879],
-                          [139.69888939,  35.53429679],
-                          [139.69889239,  35.53467279],
-                          [139.69908039,  35.53467279],
-                          [139.69907902,  35.53417],
-                          [139.69907502,  35.533637],
-                          [139.699075,  35.533637],
-                          [139.699075,  35.533637],
-                          [139.698311,  35.533642],
-                          [139.698311,  35.533796]])
+        truth = np.array([(139.6983434498457, 35.53388784954066),
+                          (139.698343456533, 35.533887872006716),
+                          (139.6984440527688, 35.53417277873741),
+                          (139.69844406145836, 35.534172800060766),
+                          (139.69844408446318, 35.534172801043766),
+                          (139.69881597541797, 35.534171000119045),
+                          (139.69881599883948, 35.53417099885312),
+                          (139.6988159998281, 35.53417097541834),
+                          (139.69881400017155, 35.533885024533475),
+                          (139.6988139988598, 35.53388500115515),
+                          (139.69881397546646, 35.53388500014851),
+                          (139.69834347324914, 35.53388784822488),
+                          (139.6983434498457, 35.53388784954066)])
         assert np.allclose(area, truth)
 
         pybdshadow.show_bdshadow(buildings=buildings,