-
Notifications
You must be signed in to change notification settings - Fork 0
/
s1_makenoisegrid.py
195 lines (153 loc) · 5.96 KB
/
s1_makenoisegrid.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
# INPUT
grid_dx_in_m = 10000
latmin = 35.5
latmax = 39.5
lonmin = -19.0
lonmax = -14.0
grid_filename = 'sourcegrid.npy'
surfel_filename = "surfel_temp.npy"
# END INPUT
import numpy as np
import os
import io
import pprint
from math import pi, cos, sqrt, sin
def wgs84():
# semi-major axis, in m
a = 6378137.0
# semi-minor axis, in m
b = 6356752.314245
# inverse flattening f
f = a / (a - b)
# squared eccentricity e
e_2 = (a ** 2 - b ** 2) / a ** 2
return(a, b, e_2, f)
def len_deg_lon(lat):
(a, b, e_2, f) = wgs84()
# This is the length of one degree of longitude
# approx. after WGS84, at latitude lat
# in m
lat = pi / 180 * lat
dlon = (pi * a * cos(lat)) / 180 * sqrt((1 - e_2 * sin(lat) ** 2))
return round(dlon, 5)
def len_deg_lat(lat):
# This is the length of one degree of latitude
# approx. after WGS84, between lat-0.5deg and lat+0.5 deg
# in m
lat = pi / 180 * lat
dlat = 111132.954 - 559.822 * cos(2 * lat) + 1.175 * cos(4 * lat)
return round(dlat, 5)
def points_on_ell(dx, xmin=-180., xmax=180., ymin=-89.999, ymax=89.999):
"""
Calculate an approximately equally spaced grid on an
elliptical Earth's surface.
:type dx: float
:param dx: spacing in latitudinal and longitudinal direction in meter
:returns: np.array(latitude, longitude) of grid points,
where -180<=lon<180 and -90 <= lat < 90
"""
if xmax <= xmin or ymax <= ymin:
msg = 'Lower bounds must be lower than upper bounds.'
raise ValueError(msg)
assert xmax <= 180., 'Longitude must be within -180 -- 180 degrees.'
assert xmin >= -180., 'Longitude must be within -180 -- 180 degrees.'
gridx = []
gridy = []
if ymin == -90.:
ymin = -89.999
warn("Resetting lat_min to -89.999 degree")
if ymax == 90.:
ymax = 89.999
warn("Resetting lat_max to 89.999 degree")
lat = ymin
# do not start or end at pole because 1 deg of longitude is 0 m there
while lat <= ymax:
d_lon = dx / len_deg_lon(lat)
# the starting point of each longitudinal circle is randomized
lon = xmin
while lon <= xmax:
gridx.append(lon)
gridy.append(lat)
lon += d_lon
d_lat = dx / len_deg_lat(lat)
lat += d_lat
return list((gridx, gridy))
def get_approx_surface_elements(lon, lat, r=6.378100e6):
if len(lon) != len(lat):
raise ValueError('Grid x and y must have same length.')
# surfel
surfel = np.zeros(lon.shape)
colat = 90. - lat
# find neighbours
for i in range(len(lon)):
# finding the relevant neighbours is very specific to how
# the grid is set up here (in rings of constant colatitude)!
# get the nearest longitude along the current colatitude
current_colat = colat[i]
if current_colat in [0., 180.]:
# surface area will be 0 at poles.
continue
colat_idx = np.where(colat == current_colat)
lon_idx_1 = np.argsort(np.abs(lon[colat_idx] - lon[i]))[1]
lon_idx_2 = np.argsort(np.abs(lon[colat_idx] - lon[i]))[2]
closest_lon_1 = lon[colat_idx][lon_idx_1]
closest_lon_2 = lon[colat_idx][lon_idx_2]
if closest_lon_1 > lon[i] and closest_lon_2 > lon[i]:
d_lon = np.abs(min(closest_lon_2, closest_lon_1) - lon[i])
elif closest_lon_1 < lon[i] and closest_lon_2 < lon[i]:
d_lon = np.abs(max(closest_lon_2, closest_lon_1) - lon[i])
else:
if closest_lon_1 != lon[i] and closest_lon_2 != lon[i]:
d_lon = np.abs(closest_lon_2 - closest_lon_1) * 0.5
else:
d_lon = np.max(np.abs(closest_lon_2 - lon[i]),
np.abs(closest_lon_1 - lon[i]))
colats = np.array(list(set(colat.copy())))
colat_idx_1 = np.argsort(np.abs(colats - current_colat))[1]
closest_colat_1 = colats[colat_idx_1]
colat_idx_2 = np.argsort(np.abs(colats - current_colat))[2]
closest_colat_2 = colats[colat_idx_2]
if (closest_colat_2 > current_colat
and closest_colat_1 > current_colat):
d_colat = np.abs(min(closest_colat_1,
closest_colat_2) - current_colat)
elif (closest_colat_2 < current_colat and
closest_colat_1 < current_colat):
d_colat = np.abs(max(closest_colat_1,
closest_colat_2) - current_colat)
else:
if (closest_colat_2 != current_colat
and closest_colat_1 != current_colat):
d_colat = 0.5 * np.abs(closest_colat_2 - closest_colat_1)
else:
d_colat = np.max(np.abs(closest_colat_2 - current_colat),
np.abs(closest_colat_1 - current_colat))
surfel[i] = (np.deg2rad(d_lon) *
np.deg2rad(d_colat) *
sin(np.deg2rad(colat[i])) * r ** 2)
return(surfel)
def create_sourcegrid(grid_dx_in_m, latmin,
latmax, lonmin, lonmax,
grid_filename):
grid = points_on_ell(grid_dx_in_m,
xmin=lonmin,
xmax=lonmax,
ymin=latmin,
ymax=latmax)
sources = np.zeros((2, len(grid[0])))
sources[0:2, :] = grid
print('Number of gridpoints: ', sources.shape[-1])
return sources
def run_s1(grid_dx_in_m, latmin, latmax,
lonmin, lonmax,
grid_filename, surfel_filename):
sourcegrid = create_sourcegrid(grid_dx_in_m, latmin,
latmax, lonmin, lonmax,
grid_filename)
surf_el = get_approx_surface_elements(sourcegrid[0], sourcegrid[1])
np.save(grid_filename, sourcegrid)
np.save(surfel_filename, surf_el)
if __name__ == "__main__":
run(grid_dx_in_m, latmin, latmax,
lonmin, lonmax,
grid_filename, surfel_filename)