diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 2c3cd56..5b257ff 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -27,7 +27,7 @@ jobs: strategy: matrix: os: [ubuntu-latest] - python-version: ["3.7", "3.8", "3.9"] + python-version: ["3.8", "3.9"] steps: - uses: actions/checkout@v2 diff --git a/example/Example2-facade shadow.ipynb b/example/Example2-facade shadow.ipynb new file mode 100644 index 0000000..b2a7b55 --- /dev/null +++ b/example/Example2-facade shadow.ipynb @@ -0,0 +1,174 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Downloading Buildings: 100%|██████████| 1/1 [00:07<00:00, 7.36s/it]\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ0AAAEQCAYAAACqWiFNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABCm0lEQVR4nO2deZxkVXn3v0/tS+/dM8PswzbAgKwji4giBCVKxKgJaDSuIUYSzWsSIzGvxhgjrzEaoybGIG5xjXEhboBCFI0sg7LDALLOMFvvXft23j/uqerbNdVV91bd2rrP9/PpT92667nV9/zuc57znPOIUgqDwWBwiq/bBTAYDP2FEQ2DweAKIxoGg8EVRjQMBoMrjGgYDAZXGNEwGAyuWBWiISLvF5F7ROQuEblBRDYss19R73OXiFxnWy8i8gEReVhEHhSRt+n1oyLyLX3u20XkJAdl+YCIPC0iCe/u0GDoHLLS4jRE5Hzg9Uqp19vWDSml5vXy24AdSqm31Dg2oZQaqLH+DcAL9HlLIrJWKXVQRP4BSCil3icixwOfVEpd2KB8ZwNPAo/UupbB0OusCkujLBiaOOBWKf8I+FulVEmf76BevwO4Sa97CNgmIusAROQ12vq4S0T+TUT8er9blVL7Wrgdg6GrrArRgMVmAfB7wHuW2S0iIrtE5FYReZlt/dHAZXrbD0TkWL3+buDl+vxnAluBTSJyAnAZcK5S6lSgqK9rMPQ9K0Y0ROQ2EbkLuAZ4qc038SIApdS7lVKbgS8Bf7zMabYqpXYCrwb+SUSO1uvDQEZv+3fgWr3+amBEX/dPgF9hCcSFwBnAHXrbhcBRXt6vwdAtVoVPo2r7FuD7Sqm6TksR+RzwXaXUN0TkIeA3lVKPi4gAs0qp4ar9BXgcOBl4HbBBKXVVnfPX9J8YDL3OirE06mFrTgBcCjxUY59REQnr5QngXOABvfnbWI5QgOcDD+v9RkQkpNe/Gfip9p/8GHiliKzV+42JyFZPb8pg6BKrQjSAq0XkPhG5B3gh8HYAEdkpItfofU4AdonI3cDNwNVKqbJoXA28QkTuBT6IJRDlY+4Tkd3Ab5bPq4/7a+AGfc0bgfX6mh8SkT1ATET2iMjftPPGDQavWXHNE4PB0F5Wi6VhMBg8ItDtAnjBxMSE2rZtW7eLYTD0BHfeeeekUmpNu86/IkRj27Zt7Nq1q9vFMBh6AhF5sp3nN80Tg8HgCiMaBoPBFUY0DAaDK4xoGAwGVxjRMBgMrjCiYTAYXGFEw2AwuMKIhsFgcMWKCO4yGPqNp6ZSzKZzJLIFktkiyWxBLxf0sl6Xs75fePxaXnvOtm4XGzCiYTB0hT//xt3c/vi04/3H4iFee04bC+QC0zwxrCqUUmTy1lu8kxSKJeYzefLFEgA71g+5On4qkWtHsZrCWBqGnqZUUqR0JU+WTflcoWLOp3LFJesT2QIpbd6nqvYrfxZLit86ZQMff9VpS65VLClSuQLpXJFkrlhZTum/dN66TmVd3r69sLif7dik/p7TYvHVK87m7KPGOf6IQVe/w1Qy69lv2ipGNAyeopQiW7DeqvPpAguZPPMZ/am/79gwxPrhCF++7WmrwucOF4SkFoNUrtiWck4uZPnojQ/zjTv3kNQVPlcoteVadmZTlsVwgrE0DCuFQrFEIltgIVNgLp1nIVNgPqM/l3xfflu+WH9ip1eftYXfPm0j1/788Q7d1eFMJrKUlGLvbLqj151J5QHYvm4Qn0DJ4RxYU4kcSimsqWi7ixGNFYRSilSuyEL5za7f8vPLVP5aVkCyTW92O7OpHBMD4bZfpx6TiSzj8VDjHT1mRlsa0ZCfjaNRnp52Jlq5Yon5TIHhaLCdxXOEEY0+o1AsUVSKcMBfWfe2r/yKnz5yiIVMgaLTV1cXmU7mmBjofIW1M5PKMxLrfBlmtaUBcMRQxLFoAEwlskY0+plsocjembR2kJUdYNXOMGtbZX3eWv/iZ63nlWdsQinFv9/yGMmy0y5XJKXb85YTb+lnKlckWyjx++ds5W8vXczA4PfJkoex15lJ5hkIBwgHfGQ74EdYjmCg86b+THLRN7F2KOLq2KlkjqPaNh+Xc1a9aOSLpYrJbjffy44863PRvJ9PF/irl5xA0C+85J9/1tQ1N41GeeUZmxARPn7Toyxk3HX/TSeXOsU2jLh7+LrNdCqHiDAxEO64T8GOj86Lxmx6qaXhhqlEb/SgrErRyBdLPO9DNzOXzjflnX96OsXpW0ebvv6UrdKPxUOuRaPaqlg/HG26LN1gNmU59SYGuysahS405cq9J+BeNA71SA/KqgzuCvp95IulprvzphJZRlpoW9pN1NEm2tXVlsbGkf4SjXxRkcgWWNNlv0a2UMLXYWNjxib464b709JYlaIBrb2dJxM5YiE/QX9zT9x0laXhlpnUUtFY32fNE7D8Gt3uQZlOZpv6/VuhFUujV2I1Vq1orHP5D7MzmcgiIgxHm3vgpjywNOxJrjb0maUBll+j26IxmcgxHu9sGWZT+cr/zq1oHFzItKNIrlm1otFKH/2kNhNHYs01UWZslX60iXNkCyXS+cWm1VAkyEC4v9xTM6nud7tOLmSZGOxsGQolxYIe97JhJELARfvosUPJdhXLFatWNEbizfskDi5YotFsn3mhpJhPWw/OaJPi1e89KDPJHBODXbY0kp23NABmk5ZfI+D3sXks5vi4xyeTHQl1b4Qj0RCRi0Vkt4g8KiLvqrE9LCJf09tvE5Fttm1X6fW7ReRFet1mEblZRB4QkftF5O22/f9BRB4SkXtE5FsiMtL6bR5OM82CMoe0aLTiDC0PQGq2TT2TXNqD0m9NFCvAq8uisZBlvAvWjt0ndeRE3PFxhZLi8cnuWxsNRUNE/MAnsbKi7wBeJSI7qnZ7EzCjlDoG+Cjw//SxO4DLgROBi4F/0ecrAH+mlNoBnA1caTvnjcBJSqmTgYeBq1q7xdo00ywoM5nIUiophls4R/nBaVa8pqudoX3W7TrTEz6NbFfKYBeNbePORQPg4QMLXhfHNU4awmcCjyqlHgMQka8ClwIP2Pa5FPgbvfwN4BNijay5FPiqUioLPC4ijwJnKqV+AewDUEotiMiDwEbgAaXUDbbz3gq8stmbq0crIcT5omIunWekSUcoLHrCm7c0qrtd+6x5kso37dMI+X1EQ37iIT/RkJ9YKEAs5Nd/i8vRUGDJPvGwn2hQ7xP2MxAO8MsnZzy+s8XyxfS146FA5Xss5GfIZqFuXzfA+uEIw9Fg7b+Y9Tmkv7sVmXbgRDQ2Ak/bvu8BzlpuH6VUQUTmgHG9/taqYzfaD9RNmdOA22pc+43A12oVSkSuAK4A2LJli4PbWEorzROw/BrNOkJh0Scx1qRv5bBu1x61NAbCAQYjAYYiQYaiAQYjQYYiAU7eOMxwNMgfnHckUV3R47qi1xKAmK3iBf3eueL2z2U4ciKuxcRPLBwgFrRV+HBgcVsNIWq1fJefuYXLz3T//HaTrrrcRWQA+C/gT5VS81Xb3o3VjPlSrWOVUp8GPg2wc+dO16F9rTRPwPJrtDJ4qNzt2qx4VVsa7fBp+H2ytMKHrc+hSNCq/JXlAEPRYGXfYb08EA4QaFCB3v2S6pZuZ3ne9jXc/Ofnd7UM/YYT0dgLbLZ936TX1dpnj4gEgGFgqt6xIhLEEowvKaW+aT+ZiLweuAS4UNkDEjyk1RGOhxIZTyyN4WgQEXB7l9U+jVq9J+GAj6Go9Wa3Krl9WYuBrvC1Kn8s5O+J+RsMvYUT0bgDOFZEjsSq8JcDr67a5zrgdcAvsHwQNymllIhcB3xZRD4CbACOBW7X/o7PAA8qpT5iP5GIXAy8E3i+UirV/K3Vp5UKD5alsX2duynb7JRFI+D3MRwNOh6lGvL7iIX9h63fOBLlO1eeW6n0g5HAkuHzBoNXNBQN7aP4Y+B6wA9cq5S6X0T+FtillLoOSwC+qB2d01jCgt7v61hO0wJwpVKqKCLPBV4L3Csid+lL/ZVS6vvAJ4AwcKN+y92qlHqLd7dsEfT7GAwHKoE2bjk4n+XMI8ebvr49zuIdF20nX1TEdZs6bnPc2T/rtZcDfh+nbB5pujwGg1Mc+TR0Zf5+1br32JYzwO8sc+wHgA9UrfsZ1B6XrLttO8JIPNi0aBxqYdBaKOCrzEoN8Ps9ks/CYHBCf8Uee8xoLORq5iSwKnw8ZJn9awbDXPmCoy3Puc3jXu7Ws3ezRYOL2xo5Bw2GXmZVi8YlJ6/n1M0jxMOWpz+uu9gGwgHi+s9atvr0Y6EAocDSCv8XLzq+S6U3GLrDqhaNK553dLeLYDD0HcZONhgMrjCiYTAYXGFEw2AwuMKIhsFgcIURDYPB4AojGgaDwRVGNAyGJlFKUeqDNJhes6rjNAwGgIVMnoMLWebSVja9ufRi4uxytj1rW0Fvy1f2/dKbz+boNXEmEzlm0znmUnnGB8KceeRYt2+rbRjRMKx6vnr703zg+w82dexUMssHf/Ag9+yZq6x7xembjGgYDCuZViYXnkrkGIwsrUaFUolUbnEgZCy0sqrZyrobQ08yl8pzYCHDgi2p9lAkyAuOX9vtogEw3sLkwlPJHIPhpaOdv3PXM3znrmcAiIf83P+3F7dUvl7DiIah7Xzuf5/goz96eMm6zWNRbjn+gi6VaCmtJM6aSmQPszTsJHNFCsXSihrZvHLuxNCzrKmRFOnp6fRh85x2i1bSGEwlcktmF6/FQqa5OVt6FWNpdIFEtsBnf/Y4U8kc08kc64bCXZ9gt53UEg2Ae/bO8fztazpcmsNpJQn0dDLH8evrT/s4n8k3nUmvFzGi4YBcoUQyWyBR/ZcpLF2fKZDMFXjF6ZvYuW1573nAJ/zjjYvm+vZ1A7z7JZ24k+6wbijMxECI9cNRjhiOsGE4wvqRaM/kagkFfAxFAsw3YRFMJrMMRupbGnNpZ/O/9gurUjRu3n2QXx9MkMwWSWTzJLJFElktAJlFEUhmCyxkC67zZ566eaSuaESCftYMhivpHaeTK+uhqubkTSPs+uuLul2MukwMhJsSjVq9J9WU8/auFFaVaCSyBd77nfv5r1/uaet1sg5EZtNotCIasykri7yTdAE/efgQjyyTmu/849ZyzNoBd4U1AFa362NN5EmdS+eJherP+j6fWVkvBUeiodMKfAxrNvJrlFJXV20PA18AzsDKd3KZUuoJve0qrFyvReBtSqnrRWSz3n8doIBPK6U+pvcfw8qqtg14AvhdpVTLufPu3TPH2776q44k0M3kiw332TQa41dPzQJWYt+FbIGhGmZusaTwCRVBuf7+/Xz5tqdqnnM4GjSi0SQtZY9vEEm+6pontgTQF2GlVbxDRK5TStlzuVYSQIvI5VgJoC+rSgC9AfiRiGxnMQH0L0VkELhTRG7U53wX8GOl1NU6Q/27gL9s5SYnE1le/q8/J1/szDiBbH7R0nh8MsmND+y3mkCZgm4OFbj76bklx7zhs3dQKCnSuQLJbJF0vkgyWyBbKPHzd13ARp1BLVInl0nGZTPKsMhYCwFehQbjT+ZXm2jQ4QTQ+pjz9bk+D/wPLYqGX6RjggFLmyePHkzw999/qOExd9ZJRLyQyQNaNILL95Jnco0tHENtJlro3WgoGiuseeIkTqNWAuiNy+2jlCoA9gTQdY+tkQB6nVJqn17ej9WEOQwRuUJEdonIrkOHDtW9gWiDNqfX2JsnreR7LWPv548G61gaDppFhtq0EhWabfC7r7TmSVeDu+olgAbQeVxryrhS6tNKqZ1KqZ1r1tTv6w8HfPg6mJLUbml4IxqLD12kjmikjWg0TSvjT5INLLyV1nviRDTcJIDGgwTQB0Rkvd5nPXDQ6c0sh4jUfUN7Tbaw+BANRVvvoLI/dHWbJ3nj02iWVhyhc+kc4cDy/5fV2DypJIAWkRCWY/O6qn3KCaDBlgBar79cRMI6gXTDBNBV53od8B23N1WLaAdHGrbT0ggbS6MtTLQ80nX5//Oqa55oH0U5AfSDwNfLCaBF5KV6t88A49rR+Q6sHg+UUvcD5QTQP0QngAbOxUoAfYGI3KX/XqzPdTVwkYg8AvyG/t4yjfrSvcTuW4gG/QT9rbWN5o1Po+204tOYTOQYqhPgtRp7TzqdAHoKuNBJudzQ2ebJoqUhIgxFgky1MDjL7git59MwotE8I9EgPoFmZu+bTtYf6dpMpGkvs2oiQlvtQQkHfAxFgwxGAgxGggxFAgxGAgxFFteVP7eOx5YcOxxtVTTsjtDaxmHAJzgIKDUsg88njMXDTCayro+dSubYMhZbdvtKa56sGtE4/ojBJc6qkrJeKUpZyyUFpZKipBTxsJ+/fsmJWgQsIahO/OyGRkOnG2G3NE7ZPMJ1f3yuzkYfqHy2Uj6DxcRAqDnRSOQ4ccPQsttzhRKZfLGuldhPrBrRmEvnue3xaUf7jsdDPGvTsGfXbl00bI7QgI+tY3GyxSL5omImlefAfJZcoUSuWCRbKJErlMgXFblCieceM8FQ1BrBmdLRpks+c0VS2apPvf11z9nKyZtGWrz7/qHZbtdEttBwSr/5TN6IRr8xEnNecZ0MOHNDqz0o9jbx31x3P1+5/ek6ey/ly39wFqdvGeWU993g+rqnbx1ZXaLRQrdroEEg0Hw6z9rB3pgKoFVWjU07HHX+FrHHWXhz7da02W5puJ2k9tBClkjQ35RwPTObdn1MP9PKZDy+Bg6luRUU4LWKRMN5pckXvU2CU2v0qhvsPo24S4duefj92mVmz6rHM7MZ18f0M63EaqgGQ11XUoDXqhENN80TgFzRuyZKq82TJWNP3Foa2rG3bsi9abx3lVkarcRqNHrHrKRYjdUjGi4rbtbDkOxWRSORLVDUT2U83KSlMeS+QuybW2Wi0cpI1wYvGSMafciwS0vDS79Gq70nYAkHNOfTgOYsjf1zmYpYrQZaGunawHm+kgK8Vo9ouLU0POxB8XL8SSd9GvmiaipuoV9pxafRaNzPSgrwWkVdru4eCC8tDS9EYz5dgFGIhTtnaYDl12j22H6jkaUR9FujpctBddGQX38GOHIiXvfYldQ8WT2i4bLiejnMvNXeE2je0phO5cgXS6xrwqcBVrfr6VtGmzq234iH/Hz+jWdaQhD0V6Jty+IQrJMlbTqZ45pbHl92+0rqPVk1omH9051P+9d7zZPmfBpKWQ90s4FFqylWQ0Q456hxK+dspsB8Os/e2XQl/+x8Os98Jk+hqLjqxScsObZRGgPTPOlDRIThaJDJhLOBY142TwYjAUSsCtwsC1ltabjsPQGridLsLOX9EKuhlCKTt0LoM/kS2UKx5ifARTus2SP3zKT4v9++ryIO8zo5tdM5Sd7+G8cuEfCg30ck6FvWQl1Js3etGtEAXIqGd5aGzycMhAMt5fRs1tIASzRO2jjMSCzIbMrdG6+TlkaxpJhO5phKZplK5JhMZHnO0eN8565nmEvna/7Np63KfuqWEW5vMLZoIBzgvve9CICAz8fNu+vPLVuPp6ZTHH/E0kFqbzj3SEolRSTor/g8okHrs6UUCT3GqhINyxnqLO+Jl3EaYAmWN6LRnKUBsG4w4l40mojVUEqRLZTI5ktkCkUyeetNn86Xl4ucuGGYD/3wIfbMpCsiMZ3KHWaNff0Pz+bDN+xu7GNyYMYlsoVKBvdWQsYBnpw6XDT+8uLjK8v5YmlJtr5EtsAdT0zz7DqZ9/qF1SUaLnwL3o8/CbJnpvm3dtmR1sxkQuWo0LVDYXYvk51tOZppnvzd9x7kMz9b3ikIcO3rd3L9/fsbxi9MJ3NsGIny2KH6Yp91GMG7kCkwGg8RCvgYjgab9jU8NZUC4P3ffYD7n5mr5PJNZAssZAo1LdWj1sS56c/Ob+p6vcSqEg03DkmvR7q22oNSbhP7fEIs5CflIsfJYqzG4c7QaNDPWDzE+ECI0ViI8XiIsXiIsYHychilFLliyWoSpKxmwWwqz2y5mZDKVZaHo0G2jtfvfnTDVDLHRgeikcw6+z3sGdwnBkJNi8aT01Z5dj0xzd175hrsbeG19dotVpdo9PHw+OqRrs2IxhvO3cZvnbLeEoV4iPF42NGMZu//7gMNLYcyJ6wf4jlHjzsuWyOmEjk2DEcb7uc0DsIuEhMDYX7dQIx++7SNvPPi4xioio8J+Kzu10YZ4+14bb12i67kctXrrwUuAQ4qpU6ynetU4FNABCt941uVUrc3f4uLjLgZHu/xfJteDlqLh/1MJpwfWxaNkzY2N7FQ3EVAWalU4ti1g/zj75yMiOD3CUJ5KkJrhjQBZlN5TtwwRK6omE/neeRg7RsqN08aMZnIOuqhsvdiTDiIkl0/HGF9HdFyNU/LarE02pHLVc9I/jngE1hiY+dDwPuUUj/QM5R/iMU0jS3hZl4Lz5snHs6p0SgmoJpDLYaCj7qoGHPpPBd+5Ceuzr9z6/LBY5OJbN2p9MqUlDUfxnSDuVjtlsYaB2NNGll0bkQjs0IsDSdjTyq5XJVSOaCcy9XOpVh5V8HK5XphdS5XpdTjwKP6fCilfgrU6iNTQPkpGQaecXE/dXETSt57zZPFN+Rg2N25ypZGs7ipGM2MVSnWMQ+mtU/DCU7EzR6Z6WSsSTJb31HrxnrNF9WKGADo5JVVKx/rWcvto5QqiIg9l+utVcdW54Gt5k+B60Xkw1ii9hwHZXSEO59G7zZP3Foaiaw1P2gzMR7gTmwLJWdvfDvFOlG6U4kcG0ediYaTZpTd0nAyqjWZayAaTYyebvb/0Cv04ijXPwL+j1JqM/B/sBIxHYabBNBlXHW5etz+9HJy4WbO1Yq1MepysJ+b5gxAobT8bz2VzHHEsLMQ+HCgsVN3vsoR2ohGvTKuB0KuAL+GE8lzk8t1j9NcrnV4HfB2vfyfwDW1dlJKfRr4NMDOnTsd2Xxd7XJtcG2/T4gF/UTKUYR6ORLwVSILlVKIiGtLAyzRaLYr1O1gv+qehkbk6lgaM6kcQZ+PNYPhhsLnJMn30t6TxhU+1cjScDsQcgX4NZz8dyu5XLEq/OXAq6v2Kedf/QW2XK4ich3wZRH5CJYj9FigUU/IM8Dzgf8BLgAecXYrjXHn0/D2n3v65lG+8gdnVwQgEvRVhMFK3ejc6Ct38w2EAwxFAgxFg5WkTdayTugUtZI5DUWDbGswdLsebi0NJ298O7k6Al0sKebSeTaMRBuKRsGBv8AeTObE0kg0sDRG492bEa5bNBQN7aMo53L1A9eWc7kCu5RS12E1Ib6oc7lOYwkLer9yLtcCi7lcEZGvYPWKTIjIHuC9SqnPAH8AfExbLBngCq9udigS4CXPWl+pXINh6/N/fz3Fvtk0uWKJbKFEOleg4HA0rFOGY0HO8Sh+4coXHM3bLzwWv5NXqwcMRgKuUha6zfRWTzSgHOAV4e4GmRucdJPbLY1RB6HkjeI/3MxyD6vH0vA8l6te/6pl9v8ZVryH5wT8Pj75e6cftv6GB/Zzz96lUX0HF3p3dKfbN3mr+HzWCOEZh+NWGs2XWU0jq24qkXUW4OVgzgq7CDjR3APzmcp4lVq4cYQevSbuqJu31+lvN65HRGpUwpWUp8ILXnDcWvbPOxNSt47QTL7E87evWTIbViwUqEyEs2EkumyAV8jvYzgWZDgaZN1QmKem64/vcTuDVqGkOLCQXbbb16lP44Lj1/JPl5/qyYRM3caIBhCukVR5LtV8wuaVyEcuO9Xxvnc+OcP37t3veP+iUnz+jWfW3eeiHevYNBplOBpkOBZkJBpiOBokEvQhuj2klOKE9/yw7ojYZmbQ2juTXlY0An4fg+EAC3XiOa58wdG846LjOtacbDdGNFjO0lg5My2tBDaPxdhcJzM7WBMtHTEU4Qk9ArUWc+l8pRfKKXtmUpx55PJD2odjQRayBQYjAYajQUa05TMcDXLJyRt48bPWO75WP2BEAwjXGG6ezBXJF0uuejUM3eeI4cNFI+ATRmIhRmNWhc4XFaGA4BPhoh3rKhW8+m9If65vECfyg7efRywUWDGWRCOMaACRGs0TsN5KTrrlDL3DX7zoeNK5IiNaIEZiIeIhf03LIhL08++/v7Pla7oZ6boSMKLB8r0RsykjGs0wMRDi9c/ZRjzsJx4OMBAOEA8FGIjo5XCAgbCfgXCQeNjvaVj1GXUGvxm8wYgG9S0Ng3u2jsf5m5ee2O1iGNqEabBjmam1mEubHhSDoRojGkAkYCwNg8EppnkCHLVmgEtP3cBA2Gp3D+p2+LOanOnKYFjJGNEAnrd9Dc/bvqbbxTAY+gLTPDEYDK4womEwGFxhRMNgMLjCiIbBYHCFEQ2DweAKIxoGg8EVRjQMBoMrTJyGwdAkpZIilS9WssUvZo7Pc9qWUdYNOUu90G8Y0TD0JaWSIlsokckXyRSKZPMlMoUimXyJXKHEfDrPQjbPQqbAQqbAfDrPfKbAQqa8zvr809/Yzm/sWMv//npqsfJnCizYBCCRtc6xKAr6M1dYNnfsYCTAX734BC7buRnfCptno+cSQOttfwJcqY/5nlLqnc3eoGFl8f179/GnX7ur7gzmx60bZPeBBUfne3omRb6oeMNn7/CqiICVEe+qb97Lt3+1lw++/FkctWbA0/N3k55LAC0iL8DKAXuKUiorImtbvUlD6yhlvdkXlryNrbd1+e27kMnb3tCFyrYFvS2RLXDcukG+9ofnNF2OkViwYcqDRIP8q3amkzniIT+hgK/heZvhtsenufhjt/D2C4/liucdtSJmgnNiaVQSQAOISDkBtF00LgX+Ri9/A/hEdQJo4HGdF+VM4BdKqZ+KyLYa1/sj4Gp9DEqpg67vyuApZ/39j5hO5sh7kAum1ZHDawcbT4o072JKg8lEFhFhIh7imbn2pK3IFUr8w/W7+e+7n+FDrzyZkzeNtOU6ncKJ7NVKAF2dxHlJAmjAngC60bHVbAfOE5HbROQnIvLsWjs1k8u1myilmM/keXo6xb175vjpw4eY7ZMZz/NF5YlgwNJE1s2wZrCxc3EhWyTg0I8wlbD+B2MOUjS2ykP7F3jZJ3/OjQ8caPu12kkvOkIDwBhwNvBs4OsicpRSS11OzeRydYpSiky+RDpfJJUrkMkXSedKrB+J8O+3PLbscc/aOMwzs2ke2r/AbCrPbCrHbDrPXCrPbDpPsSpF2Tff+hxO39L+h7VVBsIBV1ng69EoC3sjhiIBwgFfw1y7o7EQhxKNk16X72ss3plpHUsK7np6hot2rOvI9dpBLyaA3gN8U4vE7SJSAiaAtpsTH/rhQ3z250+QXia93/fe9lz+7SfLi8YlJ69n31yGO5+ccXS9Xz01W/fNOxINcsrmEUfnaiduEzrXI5EpuE4hYEdEWDsU5ukGSZEGIgFHojGl95lwkKLRK56sk2KhH+jFBNDfBl4A3Cwi24EQMOnsdlpj/Uh0WcEAywLx++Qwi6FMOlckFnKeMvH9332g7vbTtozwrbee6/h87aKZLPXLUdBdpctNseiENQONRSPu8P8wmcyhlGKsg6Lx1HR/i0ZDn4b2UZQTQD8IfL2cAFpEXqp3+wwwrh2d7wDepY+9HygngP4hhyeA/gVwnIjsEZE36XNdCxwlIvcBXwVeV900aQd3PT3LLQ/XN2YWMkWG66ThS7kUjUa4TSHYLrwUDWjdr7HWgV+jVta8WuQKJRLZQkd8GmVWg6XR6QTQOeA1TsrVKolsgX2zaX7y8CE++IOHlrUg7PuPxILLtu9T+SJHNEis44b5FiuXV3jZPAFIZguscdALshxrhxof6/c579qcTuaY6JBPA6wepLlUnmGXOW97hV50hHaETL7ISe+93tUxiWy+bsLfVLZAdAVaGgMeWxpu4ihq4aTb1Y3HZDKR62jzBODJ6SQnx0Y6ek2vWLWiEQn6mRgIMZlw3iuwkCkwElv+4UrlisRaaKtXUw6TbqX97wUDYW/fiK13uzYWjUZWo52pRJYJl5ZPwCdEgn4iQR/hgJ9w0EckYH2PBP2EAz69vdY+/rrN3F5n1YoGwMaR6GGiEQn6iAb9RIN+IiF/ZTka8jMxEK5raaTz3vo0wKpg3RYNr30ayZYtjcZNQDfRnVPJHGcdNc4/XXbq0goe9BNZsrwoCIGqyM5C0eqiT+eKuqve+svkiwyEAz3RC+YVq1o0rnrxCUwlsvhEQBRKCYWSIl8okStaA5/yxRJZ/TkcDdZthyazBaIephgEmM/kW2r/e4HXotFq88TJ75Gq0wsGEPQJAxErReQzs2nu2ztHJOgnkS2wfy5jhcHrkPhktrBk0NqnXnMGn/rJr/nxgwcrQpErLi9S5x07wRffdJbr++xVVrVofPHWJ/nePfsc7/+m5x7JSHT55km2UPLc0ugFv4bXjtCFVi0N7QgdDAcY1omeR2MhhqM66XM0xMRAiJjOHzsYCTAYCTIYCfDh63fzowcPkC8qZlJ5ZlJ5Pn7To3z8pkcdX//gQpZ8scT+eWdh5+lcfQHrN1a1aLhtV86l82wZi9XdJ+51BeuBHpR29J64JVcosZCxhrfPp3N87g3PJp0rMp/JM58uLG7L5Nk/t8C8Hvo+n86zbSLGl958NgDRoL/lkPipRJZxF4nBU0Y0Vg5DEXeiMZvKM9KgmywW8nYU43ymByyNquZJyO8jHPARLrf/Az7Cuq1fvRypLPv1/j5O3+Iss/vDBxZ4zTW3MZ/Jk8k3PwLVbv0NeeCAnE7mGHfR25JqMXS+11jVouHW0phP5+v2ngBEAl43T7r/wO3cOsZd77moIhCdmlRmKBLk4ELjUPBG2N/0Qx74Z6aSOTaORJu6/kpgVYvGUNTd7c+mc0t6T8IBH8PRIEPRIEORAMPRoOuuu0b0gqURCvgIBTo/sG7coyhN+yA5LyyNqUSOkzc5z/NrfBoriG3jcc47doKBsOVFH9B/8XCASECYSedRygoUUlifP3n4EK9/zjaKJUW2sNi1lswWmEzk+Itv3ONpGbvtCFVKkc4Xl051V+5VyJXXFUlk8ySzRRbsvQ1Za/mbb32O66YgQNDvYyweanmE7RJLw5PmSZYJNz6NfLGlQXq9xqoWjXOPmeDcYyZqbnvNNbfxs0fdj5MLeDwxU6uWRrZQrDgEF7SjcHHOzPxh2+bTtn0yeZLZAi7ipGqiWpgQa2KgddEod50H/b6mxKuaqWTOlRVU9GCQXi+xqkWjHqUmx8gVShD0i2eT1jTyaXzkxoc5tJDVFf5wYWg070QnUDT/W6wZDPPwgUTLZUjligxHfa6bpLWYSuQYs/m2Aj5hMKK7d8NBBiIBhnQ370A44HmcS7dZWXfTI8TDAWZT3jQrGlkaX7n9KQ554CxsJ8lckZH6PdXLssZFM6AeqVzB8j95YGlMJrIE/D5uf/eFDIaDRIK+FdP0cIIRjTYQC/m9E40GPo2BcKDnRWPGZW+DnUbRnwGfVBzRQ1oUhqIB/bm4vhxr4sWYj0MJK7jLSTj7SsSIRgOGIgHG4iGmkjnHgVZetl0bDY+Ph3u/nTzTwlyoF5+0ni1jsWUEwf1b3gtLQynYP5dhc4NAv5WKEY1l+NRrzyAW9FcGJr3li3fyw/v3Ozo24qE3tJGlEfN4rEs7aMWRecbWUc7Y6iwYzAluhvkvjmRdHL0a1csrLfbCDb3/xHWJVt5IXua2aOTT8DrEux141VTzAr9P+OSrTyfoXxSEqE0U7AKxEnKUtIPef+J6hJefvpGzjhqzAp38Vqi09ekj7PcR0qHSoYCPT9z0CHfvmfPkupl8iWyhSHiZSFOvx7q0A69mMveKl5y8vttF6Gt6/4nrEV544hGO9/V6GreFTIHwQG3RGFjhPo1+5NM//TVP2OYBHYoE2DubIZW1AuKS2SLJXIGhSJBvX9n9iaPdYkTDAzI6YnIhU2D9cMRzP8N8Or9sBGI/+DRmeqh50glufOAAdzyxmMbi5E3D3FPD8my2R6nb9GQCaL39z4APA2uUUh1JYbAcH//xIzw+laxkG1+SpzRTWDIBy5fefJb3c2rU6UHph+bJTI81Twyt0XMJoPU1NwMvBJ5q5ea84jt3P8OjB51FJe6byxBvg6WxHKZ5Yug0vZgAGuCjwDuB77i7ne5zYD7TcM4Nt9TrQTGWRu/xmrO38orTNzEaDzEWDzEQDlArlCTgIs1CL+HkiauVxLl6wsMlCaBFxJ4A+taqY+smgBaRS4G9Sqm76wXtiMgVwBUAW7ZscXAbreET9DgCa9q4gbA11qA8tsAaIRvk9C2jHHA4DZxT6o0/8dqqaQf95NMoFEuV5md5DE8iU2Ahm9fN08X11rbF5YtPOoL/c9F2Lj21UY7z/qannjgRiQF/hdU0qUs7E0BX890/eS7hgPPIwxscBoHVwydwwfFWkuB6Iyr7wdJI54sdScVQHtFb9jWVpwBcsFfympV/0U/VStCWV/N/9Dq9lgD6aOBIoGxlbAJ+KSJnKqVar4kOUEqR0zOQZ3WMRLZgzUxurbO+W3/Fmuu9mDgnHPBzzet2NtyvH8LIwfJrrB/2vrfgr799Lz+8bz/zmYKrtAXtYPf+hbrblbKGyJfnG8kWSmxfN9ih0nlHTyWAVkrdC6wtfxeRJ4Cd7e49ueILu/jFY1MVAWiVTR3sSuuHiFCwArzcisbe2TT//KNHSFXyiRRI2/KJvPqsLRRLylXCq3Yylczx5/95d6ULPmlLe1AWCvuUCYPhAPe+70VdLHFzNHzitI+inADaD1xbTgAN7FJKXYeVAPqL2tE5jSUs6P3KCaALHJ4A+nxgQkT2AO9VSn3G8zt0QLZQ8nTW76lk50ad9kPzBJoLJfeL8LVdTy+7fSFTYLyDOVid8I079zjedyFriaCXqTw7Qc8lgK7aZ5uT8rWK13EV6XyJgXCg5aRATugHRyg0F0q+ZjBMwGclsFqOiT73I0wmsn03WrY/nrg2E22Dg24sHuqMaPSJT2O2iVgNv0942WkbrXSXOjVmNOQnFgwQC/nZuW2UPTPpNpS2cxxcMKLRl7TDPPRiqnwn9I+l0Zxz+MO/c0rd7ekG6Rd7nV6fQKkW/fHEtRmvmyfQ3EQ8kaCvEgcy4bCt7vMJsZC/5+d3aFdUqJtZwXuRyYQRjb6kHc2TiYEwLz1lQyWX6JAWg+lkjky+SMDnw+eznH0iICIUS4qc7t7NFUtcc8tjvPm8oxpeKxYK9JxoTAyE2DgaY/NolE2jMZ67zKzvrV+nv0XDWBp9iteZ3gGOXTfAn73wuMPWX/GFXdzwwAFH53juMROORGMg7GeyiQm7Q35fJaq1HOE6WM4Bo5erI1/L+WEiQT9KwUs+fgv2idv/8PlH8fYLj/Vs9O2dT87wk90HK99/99mb2TQa48PX7+bGBw7wfy85AZ/QcpqFbnHIWBr9iZPmScjvIxryE9fOuHg4QDSoP/X6WMhy0MXDAc46cqzmeS4/czPPP24NocrEPT49sY81gU/AL2TzVrLjRsmmy5x7zATHH5E7PKzdLgjhpVPsx8N+BKnED9jjCuYriZatiMl985nDcqW8/9KTOOfocYajwcO6U70crv/IgQX+2ZbR/bQto2wajZHMFdh9YIF9cxnG4uGeMfPDAR/nHD2+5HePVwmu/f/Sj5aSEQ3gedvX8IlXn0ZcV/pYKEAs7F9cDi2d+q1QLJGuBBxZwUaV77kiqXyRxw4luW/vHOl8iXSucPh++cUgpVRucZ1fpPL2+buXncTx64calv/vXnYS371n3xIBOLSQ5YmppJX9LKOzn9mCjBLZ2hGUkYCPjIMAt8cnk5xz9Dhj8dAS0ZjyONCqevbw2XRuyfq5dJ6JgVDPiMaGkSife8OZ3S5GWzGiASQyBe58cqZuZbYiEQtk8qUl82d4zbaxGIf0stPYBhHhPd+5z5OBYSOxIPvnG1fAvbPWzFRjsRCPkaysn/K48lanUZzT9ziyRDTCQP0QboN3GNEAnppO8dmfP9HtYgAQs0V4uqmAG0ainojGQCQIDkTjmVlrJO9oPEQo4GM0FmQkGuIIj8eXHG5paNHQGc7m0vlVM1CsVzCiQW8FSEWCi82gKRdRlBtGotz/zHzL14857Enaq4Oq/vX3Tq+keWgH1bPCz2nRWD8c4di1A0RD/p6YlyIa9DMQCfTtFH5uMKJBbw36slcAN6HXXj2sIYc5W/bOWqLRTsEAK5T8/S87ifF4iPF4iK3jcQDOOmqcG9/xfAA+efOj9U6xLEG/2ByTQQbCfv09qB2Vfmt9pHp5qUMzHvK3/XfoJXqntnSRXhr0ZZ+yw41T0SvR8PuczRmyfz5DoVhqe2WJhvy89uytdfc5ccMQv3nSEYf3GNXoPbIvL5cWwlCf3qktXaSXLA17tnq3zZNOUiwpDi5kO37dWpx/3FrOP25t4x0NnrB6bKo69JKlYe8GnUnlKDmMWtow4k0y4kLJec9QuYliWF0Y0aC3HKH2AVjFkqo4/hqxcdSbN34m70I0+nyEqaE5jGhgTa0X6hFHVvVweqdNlIl42JN7SLsYw2IsjdVJ79jlXSYe9pNLdXeOSVgMXirjtAfF5xPWj0R40pYOsBnKM5iFAj6GIouzrw/qEPTBSIChqPV52paRlq5l6E+MaGji4YCnU+2HAz5iISsjeXnyGGs5QDToIxYKVDKWx/R267uPgE8qM59vctHs2DAcrYhGJOhjKBJkKBpkSFd063vAtt76PhhZ3GcwbH22e+ZwQ/9iREPz7G1jbB3PLKnYi8tWxY6EbJU8aBeCxc9YyE8k4MfnsOvSSz72qlPxiTBouhMNbcSIhuajl53a7SK0zNpBb3pQDIZ6OPKcicjFIrJbRB4VkXfV2B4Wka/p7bfZ0y2KyFV6/W4ReZFt/bUiclBE7qs61z+IyEMico+IfEtERpq/PYPB4DUNRcOWAPo3gR3Aq3RiZzuVBNBYeVj/nz7WngD6YuBf9PnASgB9cY1L3gicpJQ6GXgYuMrlPRkMPU+ppEhkC+yfy/DowQSPHOifUbo9lwBaKXWD7eutWMmXDIaWyeSL/OjBA8wkc0wlc0yXPxM5ZlI5zjl6nPf+1onLHq+UIpUrWpMVlecmySzOT1L5yxSW7pMtVOY1sSY7KpLMFZbMeHbkRJyb//z89v8IHtBzCaCreCPwtVobOp0A2tD/zKRy/PGXf7Xs9h01Jjx65b/+L3tn05YQ5Aptm1bwwHwGpZTjfMHdpDcimmogIu/Gysr2pVrblVKfVkrtVErtXLNmTWcLZ+hLGg3Gy9aYXOmxyST75jIsZNsnGACpXLEjeXK8wIlouEkATYsJoNHneD1wCfB7Sqk+nTLW0GuMx8N1haOZhE5ecrBPZibvqQTQYPXUAO8Enq+Uai280dAx9s2leWIyRTJrmfHltv2GkSi/dcqGbhcPsCyNdYNhnpnL1Nw+02RCJ684MJ/h6DUDXS2DE3oxAfQngDBwo27f3aqUeouXN91ulFKk80Xm0nnm0wWKJcWODY0nCO5nvnzbU3z8psMnw3nRiet6RjQA1o9ElxWNsqXx4L553vS5O4D2JXmqxUEH0yz2Aj2XAFp323aNeh7yZK68XGQ+k9eioD/LU/zraf7zxcVW1Ukbh/jun5zXxbvqHo8dSvLRGx8mUyiSzZfI5Iv6r0SmYC0/b/sa3np+Z/7t64eXD4ArDyMYjgaXFZZ2cmC+89dsBhMRCvz4wQO8+1v3VUzrdjq8coUS6ZzV5ZbKWd1vKT3TeTJXJJW1PtNV31O5grWf7q57x0XbufCEde0rqKao4wmS2aXdiElbN2IyV+TWx6ZqHv/IwQQf+/Ejda8x3sHcH/VEI60FrV25SOIh/9IBgJVla+zPSRuH23JdrzGigTWic38bVf6hfQuc8r4bSOUKSyyQVnjsUJILT2j9PP9x65Pc/vh0JZ7Abk0ls4WOJFiedzhniBesbzBb+kwqx/rhKKOx4JIBjOUKXx7gN1g1Anioxmjg8vahiDW3qNOpFHsdIxq0Px9owcVkOk7xypS9Z88s1939jCfnapbycPxO0GiGs5lknvXDUf7zLc8hHPB5VuGVUiSzOntd2vqcS+X19zw7t40ZS6Of6MfUeF5ZRmMOs9O3k/lMey0NpRQL2QIzDuYmKTtDj1k7QKmkWMgUeHo6xXQqx4yOIp1J5Vg3FCEeCvCThw9ZPhrd7MxqP01a+21++7QNfOZnTzCfzlOo0+5958XHGdHoJ0ZjQUSgnyJCvLI0xuLBxju1mfm0c0tDKUUyV2RGV97pymfeqtRVlXs6mWc2latU2AuOrz8BcblJ8tXbn+Ld376P4jIV/bnHTHDO0eN88dYn654vnSs6mkipOh9uL2NEAyt3x1gs5Gr2726z0i2Nm3cf5KYHDx4mAjPJfEtpMcMBH6dsHqlMkhTTkx9Fgj78YiXkBitz3HKCAVZO2fF448xuThM5ObGCegUjGprxgf4SjQPzWU/GKvSCpZErWF2x9tnCnppKNXyLN8MD++Z59raxilXyxGSSmVS+4nP6+aOT/MaOdRwxVN/3MZvKO2rWOk3+5uWsce3GiIZmzMFbo5fIFUrMpvKMtlju0Vhv3PdCprBENI5Z257IyCenUnXnUS2/OOp1zYI1l6uzHLLORL3bIexu6NkBa52m30QDvGmijPdA8wQOb6Ic2ybRaMRcOk++WGJ8IEygTo/JQrbASKyxlVZy6CjrZORpqxhLQzPSI29cN/zHrU+ydzbNM7Np/vtPntvUvKCjPdA8gcNjNdYMhhmMBDraHVtmJpVj7WCEdUORumkagnVSRhw1EWftUJj1wxFHTnbjCO1DxvpQNB7cN88vn5oFrLyvzaRIHAgHCPl9LTkXvaBaHESEY9cOVO6vXUSDfo5dN8BILMRINMhILIhP+4mOGK4vGtlCiXjIT7JGrpg3n3cUrz7Lmufl77//UMM4ndl0vm/m0zCioWnVN9AN7J75Q03mVRURRuNBDnR4sJRPrLQR5YTMvhqV5Zg2ikbI72MsHuIFx6/hgy8/ueY+Tpyh4wNhktP1B2OPxUMNRaNYUsxnCgxHe8Pyq4cRDc2og/Zpr2GvZ5OJ5iv9WDzsSDQCPllS0eNhPwORIANhP/HQYlZ2+z5Lvkf0MWErLUSjt+qxawcd30Ms5GcsHmI8HmIsHmIsHmZ8oLxsrR+1bR8IBxpe/4hGztB0jomBEE81EI1Gvg8RGIoEWcjkjWj0E/1oadidbIdamMDljeduYzaVJ66FYDASIB6yKvtgZLHShwO+jprPJ20c5vQtI5YAxEOMDdhFIcR4PFxZ147kTo16UJbrdh2JWWNOyrzunG1ccnKukpBq2J64KhpkIBToSp6cZjGioelHn4Y9LLkV0fidnZsb7+QApRRZHU6dyhdJ5xZDqtM5azSv/Xv587Jnb2bzWOyw851z9DjffOu5npStGY47YpDnHD1ONGglyopVJcbavm6QtYMRLjllA+sGwxwxbDlPqwXsZae5mRa39zGioemVeAU3ZG0Z3g+10Dxxyxd/8QRf37XnsMrf7IjYHRuGaopGtznv2DWcd6yZf7YaIxqaXul6dEMqt9jj0Iql4ZY9M2nu3Tvn2fmebuATMPQWRjQ0A+EAQb94Nt9FJ7DPXt1J0VgOEasLMxrUyazLpnywnAfXVzHty8mvo0E/z9rUH6M7DRZGNDQiwmgs1DczQgNLuvE62Tx57TlbeempG5Ykv44E/R13lPYr5dnQFjJ5FjIFFjIFEllreT5T4LnHTHDkRLzbxVwWIxo2ekk0RDis27Lcq1Hu2iz3aCiojM7sBJtGY2wa7djleop8saQrermS50noil8RgSpBWCoOhYb5TT70ipP7XzR0WoGPYc1Gfo1S6uqq7WHgC8AZWPlOLlNKPaG3XYWV67UIvE0pdb1efy1WbpODSqmTbOcaw8qqtg14AvhdpdRM03foglb9GvUq+mL3pZ+BsBXbMFBDBMrxDE7iGAzes38uw0du3H1YhZ/X1kAm3/7I2U5ajc3QUDRsCaAvwkqreIeIXKeUsudyrSSAFpHLsRJAX1aVAHoD8CMR2a7TGHwOK13BF6ou+S7gx0qpq3WG+ncBf9nKTTplx/phMvnS4ZXeVPRVg98nfH3Xnq6WoRf8U/XouQTQ+pjz9fLngf+hQ6Lxnt/a0YnLGHqY8Xio62Nxel00nDSEayWAro5WWZIAGrAngG50bDXrlFL79PJ+oOY8/SJyhYjsEpFdhw4dcnAbBkNjfD5pGD7ebnq9edLT82noPK41+0BNAmhDu2gUPt5uJleApdHpBNAHRGS9Ptd64KCDMhoMnrGxidHCXrISmieVBNAiEsJybF5XtU85ATTYEkDr9ZeLSFgnkG6YALrqXK8DvuOgjAaDZ6xvkBulnQT9QjjoI1tof5KqZunFBNBXA18XkTcBTwK/6+kdGwwNaJSFrR5+n9QezRoJ2tZZ2+3bhqPW9n4IkBPVT8k+lmHnzp1q165d3S6GYYVwyyOHuOqb9zqs9OUKb+0TC3W/211E7lRK7WzX+U1EqMFQxXnHruFnf3lBt4vRs/R074nBYOg9jGgYDAZXGNEwGAyuMKJhMBhcYUTDYDC4woiGwWBwhRENg8HgCiMaBoPBFUY0DAaDK1ZEGLmIHMIap+KGCWCyDcVpFlOe+pjy1Mdenq1KqbbNF7EiRKMZRGRXO+Pz3WLKUx9Tnvp0sjymeWIwGFxhRMNgMLhiNYvGp7tdgCpMeepjylOfjpVn1fo0DAZDc6xmS8NgMDSBEQ2DweAOpVTf/AEXA7uBR4F31dgexkrp+ChwG7DNtu0qvX438CK9LoI10fHdwP3A+2z7f0nvex9wLRDU68/HyutyF/BrrL7xTpTnc8Dj+rp3Aafq9QL8sz7XY1ipLDtRnltsZXkG+Ha7fx/bNj/wK+C7tnVH6nM8qs8ZqnGNh3SZOlGejj8/Dcrj5Pm5Bzi9YT3sthC4EAy//pGPAkL6Qd5Rtc9bgU/p5cuBr+nlHXr/sH64fq3PJ8CA3ieo/zFn6+8v1tsF+ArwR7Z/+ne7UJ7PAa+s8bu8GPiBPn6PfiDaXp6q8/4X8Pvt/n1sx70D+DJLK8XXgcv18qds/6+36u9+4IAuWyfK0/Hnp0F5Gj0/ApwN3NaoLvZT86SSHlIplQPK6SHtXIqVyhGs9JAXVqeHVEo9jqWqZyqLhN4/qP8UgFLq+3q7wnrbbupmeepwKVY+3DOxrIEoVna7jpRHRIaAC4Bvt/v30dfbBLwEuMZWBtFl+IZe9XngZVXXOBOrop0N5NtZHujO81OvPHW4FPiCLuqtwEg579By9JNotCU9pIj4ReQurKRMNyqlbrOfUESCwGuBH9pWn4P19tghIid2sDwfEJF7ROSjIhKuukb5s3yujvw+WJXzx0qp+U78PsA/Ae8E7MlWx4FZfY7q/e2/z1O2a7SzPBU6/fw0KE+956fWuWrST6LRFpRSRaXUqVhvgjNF5KSqXf4F+KlS6hb9/ZfAVuDPsd7s3+5Qea4CjgeeDYzRoaTYDn6fV2FVgDJt+31E5BLgoFLqTq/O2QoOy9Ox56dBeTx7fvpJNNqaHlIpNQvcjOVsRZ/jvcAarDZieb95bbLvxWoHBkVkot3lUUrt0yZkFvgs2hy1nav8WT5XJ36fCV2O73Xo9zkXeKmIPIFlzl8gIv+hjxnR56guv/332WK7RjvLU/59Ov38LFseB89P9bmWp5HTo1f+sHK0PIbl+Ck7jk6s2udKljqOvq6XT2Sp4+gxLEfUGmBE7xPF6hG4RH9/M/C/QLTqGkfof3ZA/7h7O1Se9fpTsEzQq/X3l2A5ssrluasT5dHr3gJ8vlO/T9Wx57PU0fefLHWEvtV+DV2eA1gC14nydPz5aVCeRs9P2RF6e8O62G0xcCkcLwYexvIWv1uv+1vgpXo5oh+eR7GcT0fZjn23Pm438Jt63clYXVP3YHWNvce2f0Hvf5f+e49e/8dYZuXdWF14T3WoPDcB9+r1/8Fir4YAn9TnegJrioC2l0dv/x/g4qp1bft9GlSKo/Q5HtXnDNe4xm6sbsdOlKfjz0+D8jh5fu4FdjaqhyaM3GAwuKKffBoGg6EHMKJhMBhcYUTDYDC4woiGwWBwhRENg8EFIvI7InK/iJREpOacnCISEZHbReRuve/7bNsuEJFfish9IvL5cmyJiAyLyH/bjnlDg3JsFpGbReQBvf/bvb3T5TGiYTAsg4icLyKfq1p9H/By4Kd1Ds0CFyilTgFOBS4WkbNFxIc1luRypdRJWN3jr9PHXAk8oI85H/hHEQnVuUYB+DOl1A6s+IorRWSHm/trFiMaBoMLlFIPKqV2N9hHqdoD/caBnFLqYb3tRuAV5cOAQT0gbQCYxhIGROQvROQOPW7kffoa+5RSv9TLC8CDNBgz4hVGNAyGNrDMQL9JIGBr1rySxRDuTwAnYM1Nci/wdqVUSUReCByLFfZ9KnCGiDyv6lrbgNOwpi5oO4HGuxgMqwsRuQ0rRHsAGNOVH+AvlVLXOzmHUqoInCoiI8C3ROQkpdR9InI5UB5legNQ1Ie8CCty9ALgaOBGEbkFeKH++5XebwBLRH6qyzqANZ/Jn6qlI43bhhENg6EKpdRZYPk0gNcrpV7fwrlmRaQ80O8+pdQvgPP0+V8IbNe7vgFrPIgCHhWRx7FGpQrwQaXUv1WfWw+7/y/gS0qpbzZbRreY5onB4DEiskZbGIhIFLgIa5wJIrJWf4axhqd/Sh/2FHCh3rYOOA5rINr1wBu1RYGIbBSRtdr38RngQaXURzp0a4ARDYPBFSLy2yKyB2sine+JyPV6/QYR+b7ebT1ws4jcA9yB5dP4rt72FyLyINYgwP9WSt2k178feI6I3Av8GKspNKmUugFr6r5f6G3fAAaxhsG/Fmv4+13678Xtvn8weU8MBoNLjKVhMBhcYUTDYDC4woiGwWBwhRENg8HgCiMaBoPBFUY0DAaDK4xoGAwGV/x/Q9WsC89puY8AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "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": null, + "metadata": {}, + "outputs": [], + "source": [ + "buildings_gdf.crs" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "Must pass either crs or epsg.", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m/Users/yuqing/Nutstore Files/我的坚果云/python_new/2022/pybdshadow/src/Example2-facade shadow.ipynb 单元格 2\u001b[0m line \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m precision \u001b[39m=\u001b[39m \u001b[39m3600\u001b[39m\n\u001b[1;32m 2\u001b[0m date \u001b[39m=\u001b[39m \u001b[39m'\u001b[39m\u001b[39m2022-01-01\u001b[39m\u001b[39m'\u001b[39m\n\u001b[0;32m----> 3\u001b[0m wallsunshine \u001b[39m=\u001b[39m pybdshadow\u001b[39m.\u001b[39;49mcal_sunshine_facade(buildings_gdf, date, precision)\n\u001b[1;32m 6\u001b[0m floorsunshine \u001b[39m=\u001b[39m pybdshadow\u001b[39m.\u001b[39mcal_sunshine(buildings_gdf,\n\u001b[1;32m 7\u001b[0m day\u001b[39m=\u001b[39mdate,\n\u001b[1;32m 8\u001b[0m roof\u001b[39m=\u001b[39m\u001b[39mFalse\u001b[39;00m,\n\u001b[1;32m 9\u001b[0m accuracy\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mvector\u001b[39m\u001b[39m'\u001b[39m,\n\u001b[1;32m 10\u001b[0m precision\u001b[39m=\u001b[39mprecision)\n\u001b[1;32m 11\u001b[0m floorsunshine[\u001b[39m'\u001b[39m\u001b[39mheight\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m \u001b[39m0\u001b[39m\n", + "File \u001b[0;32m~/Nutstore Files/我的坚果云/python_new/2022/pybdshadow/src/pybdshadow/facade.py:572\u001b[0m, in \u001b[0;36mcal_sunshine_facade\u001b[0;34m(buildings_gdf, day, precision, padding)\u001b[0m\n\u001b[1;32m 566\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mcal_sunshine_facade\u001b[39m(buildings_gdf, day, precision\u001b[39m=\u001b[39m\u001b[39m3600\u001b[39m, padding\u001b[39m=\u001b[39m\u001b[39m1800\u001b[39m):\n\u001b[1;32m 567\u001b[0m \n\u001b[1;32m 568\u001b[0m \u001b[39m# 计算阴影重叠情况\u001b[39;00m\n\u001b[1;32m 569\u001b[0m final_shadow \u001b[39m=\u001b[39m calculate_buildings_shadow_overlap(\n\u001b[1;32m 570\u001b[0m buildings_gdf, day, precision\u001b[39m=\u001b[39mprecision, padding\u001b[39m=\u001b[39mpadding)\n\u001b[0;32m--> 572\u001b[0m final_shadow[\u001b[39m'\u001b[39m\u001b[39mbuilding_id\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m final_shadow[\u001b[39m'\u001b[39m\u001b[39mbuilding_index\u001b[39m\u001b[39m'\u001b[39m]\n\u001b[1;32m 574\u001b[0m \u001b[39m# 求最大光照时长\u001b[39;00m\n\u001b[1;32m 575\u001b[0m lon, lat \u001b[39m=\u001b[39m buildings_gdf[\u001b[39m'\u001b[39m\u001b[39mgeometry\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39miloc[\u001b[39m0\u001b[39m]\u001b[39m.\u001b[39mbounds[:\u001b[39m2\u001b[39m]\n", + "File \u001b[0;32m~/Nutstore Files/我的坚果云/python_new/2022/pybdshadow/src/pybdshadow/facade.py:511\u001b[0m, in \u001b[0;36mcalculate_buildings_shadow_overlap\u001b[0;34m(buildings_gdf, date, precision, padding)\u001b[0m\n\u001b[1;32m 508\u001b[0m shadows_gdf\u001b[39m.\u001b[39mset_crs(buildings_gdf\u001b[39m.\u001b[39mcrs, inplace\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m)\n\u001b[1;32m 509\u001b[0m shadows_gdf \u001b[39m=\u001b[39m shadows_gdf[[\u001b[39m'\u001b[39m\u001b[39mbuilding_id\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mgeometry\u001b[39m\u001b[39m'\u001b[39m]]\n\u001b[0;32m--> 511\u001b[0m overlapping \u001b[39m=\u001b[39m gpd\u001b[39m.\u001b[39msjoin(buildings_gdf_overlap,\n\u001b[1;32m 512\u001b[0m shadows_gdf, how\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mleft\u001b[39m\u001b[39m'\u001b[39m, op\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mintersects\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[1;32m 514\u001b[0m sun_vec \u001b[39m=\u001b[39m sun_light_vector(sun_azimuth, sun_altitude)\n\u001b[1;32m 515\u001b[0m walls_target[\u001b[39m'\u001b[39m\u001b[39msun_vector\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m [sun_vec] \u001b[39m*\u001b[39m \u001b[39mlen\u001b[39m(walls_target)\n", + "File \u001b[0;32m~/miniforge3/envs/py38_native/lib/python3.8/site-packages/geopandas/geodataframe.py:1279\u001b[0m, in \u001b[0;36mGeoDataFrame.set_crs\u001b[0;34m(self, crs, epsg, inplace, allow_override)\u001b[0m\n\u001b[1;32m 1277\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[1;32m 1278\u001b[0m df \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\n\u001b[0;32m-> 1279\u001b[0m df\u001b[39m.\u001b[39mgeometry \u001b[39m=\u001b[39m df\u001b[39m.\u001b[39;49mgeometry\u001b[39m.\u001b[39;49mset_crs(\n\u001b[1;32m 1280\u001b[0m crs\u001b[39m=\u001b[39;49mcrs, epsg\u001b[39m=\u001b[39;49mepsg, allow_override\u001b[39m=\u001b[39;49mallow_override, inplace\u001b[39m=\u001b[39;49m\u001b[39mTrue\u001b[39;49;00m\n\u001b[1;32m 1281\u001b[0m )\n\u001b[1;32m 1282\u001b[0m \u001b[39mreturn\u001b[39;00m df\n", + "File \u001b[0;32m~/miniforge3/envs/py38_native/lib/python3.8/site-packages/geopandas/geoseries.py:1031\u001b[0m, in \u001b[0;36mGeoSeries.set_crs\u001b[0;34m(self, crs, epsg, inplace, allow_override)\u001b[0m\n\u001b[1;32m 1029\u001b[0m crs \u001b[39m=\u001b[39m CRS\u001b[39m.\u001b[39mfrom_epsg(epsg)\n\u001b[1;32m 1030\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m-> 1031\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mMust pass either crs or epsg.\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m 1033\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m allow_override \u001b[39mand\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcrs \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mand\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcrs \u001b[39m==\u001b[39m crs:\n\u001b[1;32m 1034\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m 1035\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mThe GeoSeries already has a CRS which is not equal to the passed \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 1036\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mCRS. Specify \u001b[39m\u001b[39m'\u001b[39m\u001b[39mallow_override=True\u001b[39m\u001b[39m'\u001b[39m\u001b[39m to allow replacing the existing \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 1037\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mCRS without doing any transformation. If you actually want to \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 1038\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mtransform the geometries, use \u001b[39m\u001b[39m'\u001b[39m\u001b[39mGeoSeries.to_crs\u001b[39m\u001b[39m'\u001b[39m\u001b[39m instead.\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 1039\u001b[0m )\n", + "\u001b[0;31mValueError\u001b[0m: Must pass either crs or epsg." + ] + } + ], + "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/requirements.txt b/requirements.txt index 2402469..4e48810 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,3 +6,8 @@ suncalc keplergl scikit-opt transbigdata +mapbox_vector_tile +vt2geojson +requests +tqdm +retrying \ No newline at end of file diff --git a/setup.py b/setup.py index 6c8f66b..0a7e54c 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setuptools.setup( name="pybdshadow", - version="0.3.4", + version="0.3.5", author="Qing Yu", author_email="qingyu0815@foxmail.com", description="Python package to generate building shadow geometry", @@ -17,7 +17,7 @@ "Bug Tracker": "https://github.com/ni1o1/pybdshadow/issues", }, install_requires=[ - "numpy", "pandas", "shapely", "geopandas", "matplotlib","suncalc","keplergl","transbigdata" + "numpy", "pandas", "shapely", "geopandas", "matplotlib","suncalc","keplergl","transbigdata","mapbox_vector_tile","vt2geojson","requests","tqdm","retrying" ], classifiers=[ "Operating System :: OS Independent", @@ -26,11 +26,10 @@ "Topic :: Software Development :: Libraries :: Python Modules", "License :: OSI Approved :: BSD License", "Programming Language :: Python", - "Programming Language :: Python :: 3.7", "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", ], package_dir={'pybdshadow': 'src/pybdshadow'}, packages=['pybdshadow'], - python_requires=">=3.6", + python_requires=">=3.8", ) 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..8a2e572 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 = {'init': '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 +