-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathshadowingfunctions.py
201 lines (184 loc) · 7 KB
/
shadowingfunctions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
# -*- coding: utf-8 -*-
# Ready for python action!
import numpy as np
# import matplotlib.pylab as plt
def shadowingfunctionglobalradiation(a, azimuth, altitude, scale, dlg, forsvf):
#%This m.file calculates shadows on a DEM
#% conversion
degrees = np.pi/180.
azimuth = np.dot(azimuth, degrees)
altitude = np.dot(altitude, degrees)
#% measure the size of the image
sizex = a.shape[0]
sizey = a.shape[1]
if forsvf == 0:
barstep = np.max([sizex, sizey])
dlg.progressBar.setRange(0, barstep)
dlg.progressBar.setValue(0)
#% initialise parameters
f = a
dx = 0.
dy = 0.
dz = 0.
temp = np.zeros((sizex, sizey))
index = 1.
#% other loop parameters
amaxvalue = a.max()
pibyfour = np.pi/4.
threetimespibyfour = 3.*pibyfour
fivetimespibyfour = 5.*pibyfour
seventimespibyfour = 7.*pibyfour
sinazimuth = np.sin(azimuth)
cosazimuth = np.cos(azimuth)
tanazimuth = np.tan(azimuth)
signsinazimuth = np.sign(sinazimuth)
signcosazimuth = np.sign(cosazimuth)
dssin = np.abs((1./sinazimuth))
dscos = np.abs((1./cosazimuth))
tanaltitudebyscale = np.tan(altitude) / scale
#% main loop
while (amaxvalue >= dz and np.abs(dx) < sizex and np.abs(dy) < sizey):
if forsvf == 0:
dlg.progressBar.setValue(index)
#while np.logical_and(np.logical_and(amaxvalue >= dz, np.abs(dx) <= sizex), np.abs(dy) <= sizey):(np.logical_and(amaxvalue >= dz, np.abs(dx) <= sizex), np.abs(dy) <= sizey):
#if np.logical_or(np.logical_and(pibyfour <= azimuth, azimuth < threetimespibyfour), np.logical_and(fivetimespibyfour <= azimuth, azimuth < seventimespibyfour)):
if (pibyfour <= azimuth and azimuth < threetimespibyfour or fivetimespibyfour <= azimuth and azimuth < seventimespibyfour):
dy = signsinazimuth * index
dx = -1. * signcosazimuth * np.abs(np.round(index / tanazimuth))
ds = dssin
else:
dy = signsinazimuth * np.abs(np.round(index * tanazimuth))
dx = -1. * signcosazimuth * index
ds = dscos
#% note: dx and dy represent absolute values while ds is an incremental value
dz = ds *index * tanaltitudebyscale
temp[0:sizex, 0:sizey] = 0.
absdx = np.abs(dx)
absdy = np.abs(dy)
xc1 = (dx+absdx)/2.+1.
xc2 = sizex+(dx-absdx)/2.
yc1 = (dy+absdy)/2.+1.
yc2 = sizey+(dy-absdy)/2.
xp1 = -((dx-absdx)/2.)+1.
xp2 = sizex-(dx+absdx)/2.
yp1 = -((dy-absdy)/2.)+1.
yp2 = sizey-(dy+absdy)/2.
temp[int(xp1)-1:int(xp2), int(yp1)-1:int(yp2)] = a[int(xc1)-1:int(xc2), int(yc1)-1:int(yc2)]-dz
f = np.maximum(f, temp)
index += 1.
f = f-a
f = np.logical_not(f)
sh = np.double(f)
#sh = 1.-f
return sh
def shadowingfunction_20(a, vegdem, vegdem2, azimuth, altitude, scale, amaxvalue, bush, dlg, forsvf):
#% This function casts shadows on buildings and vegetation units
#% conversion
degrees = np.pi/180.
azimuth = np.dot(azimuth, degrees)
altitude = np.dot(altitude, degrees)
#% measure the size of the image
sizex = a.shape[0]
sizey = a.shape[1]
#% initialise parameters
if forsvf == 0:
barstep = np.max([sizex, sizey])
dlg.progressBar.setRange(0, barstep)
dlg.progressBar.setValue(0)
dx = 0.
dy = 0.
dz = 0.
temp = np.zeros((sizex, sizey))
tempvegdem = np.zeros((sizex, sizey))
tempvegdem2 = np.zeros((sizex, sizey))
sh = np.zeros((sizex, sizey))
vbshvegsh = np.zeros((sizex, sizey))
vegsh = np.zeros((sizex, sizey))
tempbush = np.zeros((sizex, sizey))
f = a
g = np.zeros((sizex, sizey))
bushplant = bush > 1.
pibyfour = np.pi/4.
threetimespibyfour = 3.*pibyfour
fivetimespibyfour = 5.*pibyfour
seventimespibyfour = 7.*pibyfour
sinazimuth = np.sin(azimuth)
cosazimuth = np.cos(azimuth)
tanazimuth = np.tan(azimuth)
signsinazimuth = np.sign(sinazimuth)
signcosazimuth = np.sign(cosazimuth)
dssin = np.abs((1./sinazimuth))
dscos = np.abs((1./cosazimuth))
tanaltitudebyscale = np.tan(altitude) / scale
index = 1
#% main loop
while (amaxvalue >= dz and np.abs(dx) < sizex and np.abs(dy) < sizey):
if forsvf == 0:
dlg.progressBar.setValue(index)
if (pibyfour <= azimuth and azimuth < threetimespibyfour or fivetimespibyfour <= azimuth and azimuth < seventimespibyfour):
dy = signsinazimuth * index
dx = -1. * signcosazimuth * np.abs(np.round(index / tanazimuth))
ds = dssin
else:
dy = signsinazimuth * np.abs(np.round(index * tanazimuth))
dx = -1. * signcosazimuth * index
ds = dscos
#% note: dx and dy represent absolute values while ds is an incremental value
dz = np.dot(np.dot(ds, index), tanaltitudebyscale)
tempvegdem[0:sizex, 0:sizey] = 0.
tempvegdem2[0:sizex, 0:sizey] = 0.
temp[0:sizex, 0:sizey] = 0.
absdx = np.abs(dx)
absdy = np.abs(dy)
xc1 = (dx+absdx)/2.+1.
xc2 = sizex+(dx-absdx)/2.
yc1 = (dy+absdy)/2.+1.
yc2 = sizey+(dy-absdy)/2.
xp1 = -((dx-absdx)/2.)+1.
xp2 = sizex-(dx+absdx)/2.
yp1 = -((dy-absdy)/2.)+1.
yp2 = sizey-(dy+absdy)/2.
tempvegdem[int(xp1)-1:xp2, int(yp1)-1:yp2] = vegdem[int(xc1)-1:xc2, int(yc1)-1:yc2]-dz
tempvegdem2[int(xp1)-1:xp2, int(yp1)-1:yp2] = vegdem2[int(xc1)-1:xc2, int(yc1)-1:yc2]-dz
temp[int(xp1)-1:xp2, int(yp1)-1:yp2] = a[int(xc1)-1:xc2, int(yc1)-1:yc2]-dz
f = np.maximum(f, temp)
sh[(f > a)] = 1.
sh[(f <= a)] = 0.
#%Moving building shadow
fabovea = tempvegdem > a
#%vegdem above DEM
gabovea = tempvegdem2 > a
#%vegdem2 above DEM
vegsh2 = fabovea-gabovea
vegsh = np.maximum(vegsh, vegsh2)
vegsh[(vegsh*sh > 0.)] = 0.
#% removing shadows 'behind' buildings
vbshvegsh = vegsh+vbshvegsh
#% vegsh at high sun altitudes
if index == 1.:
firstvegdem = tempvegdem-temp
firstvegdem[(firstvegdem <= 0.)] = 1000.
vegsh[(firstvegdem < dz)] = 1.
vegsh = vegsh*(vegdem2 > a)
vbshvegsh = np.zeros((sizex, sizey))
#% Bush shadow on bush plant
if np.logical_and(bush.max() > 0., np.max((fabovea*bush)) > 0.):
tempbush[0:sizex, 0:sizey] = 0.
tempbush[int(xp1)-1:xp2, int(yp1)-1:yp2] = bush[int(xc1)-1:xc2,int(yc1)-1:yc2]-dz
g = np.maximum(g, tempbush)
g *= bushplant
index += 1.
sh = 1.-sh
vbshvegsh[(vbshvegsh > 0.)] = 1.
vbshvegsh = vbshvegsh-vegsh
if bush.max() > 0.:
g = g-bush
g[(g > 0.)] = 1.
g[(g < 0.)] = 0.
vegsh = vegsh-bushplant+g
vegsh[(vegsh<0.)] = 0.
vegsh[(vegsh > 0.)] = 1.
vegsh = 1.-vegsh
vbshvegsh = 1.-vbshvegsh
shadowresult = {'sh': sh, 'vegsh': vegsh, 'vbshvegsh': vbshvegsh}
return shadowresult