-
Notifications
You must be signed in to change notification settings - Fork 3
/
Gravi4GW.py
263 lines (228 loc) · 11.6 KB
/
Gravi4GW.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
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
# -*- coding: utf-8 -*-
"""
Gravi4GW.py
Landon Halloran (ljsh.ca) 2021
https://github.com/lhalloran/Gravi4GW
"""
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as interp
from osgeo import gdal
# define constants (all in SI units)
G = 6.67408E-11 # gravitational constant
rho_H2O = 1000 # density of water
def quickplotterxyz(x,y,z):
"""
quickplotterxyz(x,y,z)
Simple plotting of x,y,z arrays for debugging/sanity checks.
x,y,z : 2-D arrays of x, y, and elevation
"""
fig, axs = plt.subplots(nrows=1,ncols=3,sharex=True,sharey=True,figsize=(12,6))
axs[0].imshow(x); axs[0].set_title('x, min='+str(int(min(x.flatten())))+', max='+str(int(max(x.flatten()))))
axs[1].imshow(y); axs[1].set_title('y, min='+str(int(min(y.flatten())))+', max='+str(int(max(y.flatten()))))
axs[2].imshow(z); axs[2].set_title('z, min='+str(int(min(z.flatten())))+', max='+str(int(max(z.flatten()))))
def pointmaker(stn_x, stn_y, max_r, nr, dens_azi=8):
"""
pointmaker(stn_x, stn_y, max_r, nr, dens_azi=8)
Returns x,y coordinates the element area represented by
each point. Point distribution is radial with increasing
distance from centre for each subsequent radius. Even
radial spacing for each r.
Parameters:
stn_x, stn_y : coordinates of centre point
max_r : distance of furthest points from centre point
nr : number of radii (including r=0) for point creation
dens_azi : (optional) point density at first non-zero
radius, controls the azimuthal point
density. must be multiple of 4
"""
if dens_azi%4 != 0:
print('#Gravi4GW: pointmaker requires a radial density parameter dens_azi that is a multiple of 4')
return 0
rs = np.append(np.array(0.0),np.logspace(np.log10(0.1), np.log10(max_r), num=nr-1))
log_fac = np.log10(rs[-1]/rs[-2])
rsextra = rs[-1]*10**log_fac # for calculation of elemental areas of farthest points
n = rs-rs
n = dens_azi + 4*np.floor_divide(rs,10) #number of points for given r, increasing with radius
n[0] = 1 # one point at centre
A = rs-rs # area represent by each point at given r
for i in np.arange(1,nr-1):
r_out = (rs[i+1]+rs[i])/2 # outer radius of element
r_in = (rs[i]+rs[i-1])/2 # inter radius of element
A[i] = np.pi*(r_out**2-r_in**2)/n[i]
A[0] = np.pi*((rs[0]+rs[1])/2)**2
A[-1] = np.pi*((rsextra+rs[-1])**2-(rs[-1]+rs[-2])**2)/(4*n[-1])
# make r,phi coordinates, referenced to centre point
pt_r = np.array([])
pt_phi = np.array([])
pt_A = np.array([])
for i in np.arange(nr):
for j in np.arange(n[i]):
pt_r = np.append(pt_r,rs[i])
pt_phi = np.append(pt_phi,2*np.pi*j/n[i])
pt_A = np.append(pt_A,A[i])
# convert to x,y cords:
x = stn_x+np.multiply(pt_r,np.cos(pt_phi))
y = stn_y+np.multiply(pt_r,np.sin(pt_phi))
return x,y,pt_A
def dg(xyz_stn,xyz_pt,dm,debug=False):
"""
dg(xyz_stn,xyz_pt,dm,debug=False)
Function that returns the delta gravity vector for a station point and
a grid point, given a change in mass of dm at the point
xyz_stn : local coordiates (in m) of gravity station. numpy array of size 3.
xyz_pt : local coordiates (in m) of point on GW table grid. numpy array of size 3.
dm : (near infinitessimal) change in mass of water at the point on the grid
debug : True will print internal values for each call of the function
"""
dxyz = xyz_pt-xyz_stn;
dx = dxyz[0]
dy = dxyz[1]
dz = dxyz[2]
d = np.sqrt(dx**2+dy**2+dz**2) #distance between station and point
dg1 = G*dm/d**2
theta = np.arctan2(dz,np.sqrt(dx**2+dy**2)) # angle below xy surface
phi = np.arctan2(dy,dx)
dgx = dg1*np.cos(theta)*np.sin(phi)
dgy = dg1*np.cos(theta)*np.cos(phi)
dgz = dg1*np.sin(theta)
dgxyz = np.array([dgx,dgy,dgz])
if(debug):
print(dxyz)
print([d,dg1,theta,phi])
print(dgxyz)
return dgxyz
def hillshade(array,azimuth,angle_altitude):
"""
Gravi4GW.hillshade(array,azimuth,angle_altitude)
Adapted from github.com/rveciana/introduccion-python-geoespacial/blob/master/hillshade.py
Arguments:
array: 2-D numpy array
azimuth: azimuthal angle of lightsource (degrees)
angle_altitude: altitude angle of lightsource (degrees)
"""
azimuth = 360.0 - azimuth
x, y = np.gradient(array)
slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
aspect = np.arctan2(-x, y)
azimuthrad = azimuth*np.pi/180.
altituderad = angle_altitude*np.pi/180.
shaded = np.sin(altituderad)*np.sin(slope) + np.cos(altituderad)*np.cos(slope)*np.cos((azimuthrad - np.pi/2.) - aspect)
return 255*(shaded + 1)/2
def Gravi4GW(tif_path, gravstn_x, gravstn_y, GW_depth, accept_resid=0.025, n_r=30, dens_azi=8, do_figs=True):
"""
Gravi4GW.Gravi4GW(tif_path, gravstn_x, gravstn_y, GW_depth, accept_resid=0.025, n_r=30, dens_azi=8, do_figs=True)
Arguments:
tif_path: string
Path to the geotiff file (e.g. DEM or water table elevation model).
Coordinate-system must be meter-based (i.e., not degrees).
gravstn_x, gravstn_y: array-like or scalar
The x and y coordinate(s) of the stations in the same coordinate systems as the
geotiff at tif_path. If one or both of these arrays has length > 1, a grid is formed
at all x,y pairs.
GW_depth: scalar
The estimated vertical distance between the gravity sensor location and the water
table.
Optional arguments:
accept_resid: approximate acceptable residual based on Bouger plate approximation (see paper)
n_r: number of radial distances for point definition in numerical integral (default = 30)
dens_azi: azimuthal density parameter for integration mesh creation, multiple of 4 (default = 8)
do_figs: boolean to enable creation of figures.
"""
DEM_in = gdal.Open(tif_path, gdal.GA_ReadOnly)
print('#Gravi4GW: File '+str(tif_path)+' imported. Size = '+str(DEM_in.RasterXSize)+' x '+str(DEM_in.RasterYSize))
DEM_z = np.array(np.float64(DEM_in.ReadAsArray()))
# make x,y (long, lat) matrices for input DEM
print('#Gravi4GW: Creating x,y matrices...')
DEM_geotransform = DEM_in.GetGeoTransform()
xy_inds = np.indices((DEM_in.RasterXSize, DEM_in.RasterYSize))
DEM_x = DEM_geotransform[0] + DEM_geotransform[1]*xy_inds[0] + DEM_geotransform[2]*xy_inds[1]
DEM_y = DEM_geotransform[3] + DEM_geotransform[4]*xy_inds[0] + DEM_geotransform[5]*xy_inds[1]
DEM_x = DEM_x.transpose()
DEM_y = DEM_y.transpose()
if do_figs:
quickplotterxyz(DEM_x,DEM_y,DEM_z)
# create interpolation function
print('#Gravi4GW: Creating DEM interpolation function...')
interp_spline = interp.RectBivariateSpline(DEM_x[0,:],-DEM_y[:,0],DEM_z.transpose()) # x,y must be strictly increasing, hence - sign on y
stn_array_size = [np.size(gravstn_x),np.size(gravstn_y)]
nstns = stn_array_size[0]*stn_array_size[1]
stn_x,stn_y = np.meshgrid(gravstn_x,gravstn_y)
stn_x = stn_x.flatten()
stn_y = stn_y.flatten()
stn_z = stn_x-stn_x
for i in np.arange(nstns):
stn_z[i] = interp_spline(stn_x[i],-stn_y[i]) #+ 0.5
max_r = GW_depth/accept_resid # max radius from stn for calculation (see notes)
dataproc = []
# cut out relevant part of DEM and x,y arrays:
mrgn = max_r
mmx = [np.min(gravstn_x)-mrgn, np.max(gravstn_x)+mrgn]
mmy = [np.min(gravstn_y)-mrgn, np.max(gravstn_y)+mrgn]
xinds = np.where(np.logical_and(np.min(DEM_x,axis=0) > mmx[0], np.max(DEM_x,axis=0) < mmx[1]))[0]
yinds = np.where(np.logical_and(np.min(DEM_y,axis=1) > mmy[0], np.max(DEM_y,axis=1) < mmy[1]))[0]
DEM_xC = DEM_x[yinds,:][:,xinds]
DEM_yC = DEM_y[yinds,:][:,xinds]
DEM_zC = DEM_z[yinds,:][:,xinds]
for n in np.arange(nstns):
print('Station point '+str(n+1)+' of '+str(nstns))
stn_xyz = np.array([stn_x[n],stn_y[n],stn_z[n]])
GW_d = GW_depth # this could be non-constant in future versions.
# create sampling points and calulate DEM at these points
print('#Gravi4GW: Defining sampling points and calculating interpolated DEM at these points...')
xx,yy,AA = pointmaker(stn_xyz[0],stn_xyz[1],max_r,n_r,dens_azi=dens_azi)
npts = np.size(xx)
zz = xx-xx
for i in np.arange(npts):
zz[i] = interp_spline(xx[i],-yy[i])
if n==0 and do_figs: # plot of points for integration (only do this once)
plt.figure(figsize=(10,8)); plt.scatter(xx,yy,c=zz,s=AA,alpha=0.5); plt.axis('equal'); plt.title(str(npts)+' Integration Points'); plt.colorbar()# plot points
# Calculate dg/dH by numerical integration:
dH = 0.01 # drop in water table in equivalent H2O (equivalent to delta hydraulic head / porosity). Should be small so as to estimate dG/dH.
dgxyz_sum = np.array([0,0,0])
progresspct = 0
print('#Gravi4GW: Evaluating delta g integral...')
for i in np.arange(npts):
if int(100*i/npts)-progresspct >=20: # print progress
progresspct = int(100*i/npts)
print('#Gravi4GW: Integration progress = '+str(progresspct)+'%')
dm = dH*rho_H2O*AA[i]
pt_xyz = np.array([xx[i],yy[i],zz[i]-GW_d])
dgxyz_sum = dgxyz_sum + dg(stn_xyz,pt_xyz,dm)
dgxyz_sum[2]=-dgxyz_sum[2] # switch to z-component positive convention
dg1 = np.append(dgxyz_sum,np.sqrt(dgxyz_sum[0]**2+dgxyz_sum[1]**2+dgxyz_sum[2]**2))
dgdH = dg1/dH
dgdH_uGal = dgdH*1E8 #in microGal/mH2O
dataproc.append([stn_xyz[0],stn_xyz[1],stn_xyz[2],GW_d,dgdH_uGal[0],dgdH_uGal[1],dgdH_uGal[2],dgdH_uGal[3]])
print('#Gravi4GW: beta_x = ' + str(dgdH_uGal[0]) + ' uGal/mH2O')
print('#Gravi4GW: beta_y = ' + str(dgdH_uGal[1]) + ' uGal/mH2O')
print('#Gravi4GW: beta_z = ' + str(dgdH_uGal[2]) + ' uGal/mH2O')
print('#Gravi4GW: beta = ' + str(dgdH_uGal[3]) + ' uGal/mH2O')
dataproc = np.array(dataproc) #convert data to numpy array
#%% plot the results
only1xy = np.size(gravstn_x) > 1 and np.size(gravstn_y) > 1
if do_figs and only1xy:
fig, axs = plt.subplots(nrows=1,ncols=2,sharex=True,sharey=True,figsize=(15,8))
DEM_hs = hillshade(DEM_zC,45,20)
cbobj1 = axs[0].imshow(np.flipud(DEM_hs),cmap='Greys',alpha=1, interpolation='bilinear',extent=[DEM_xC[0,0],DEM_xC[0,-1],DEM_yC[0,0],DEM_yC[-1,0]])
demobj = axs[0].contourf(DEM_xC,DEM_yC,DEM_zC,cmap='gist_earth',alpha=0.5, levels=80)
plt.gca().invert_yaxis()
axs[0].set_title('Input data (m)')
axs[0].set_aspect('equal')
axs[0].grid(c='k')
plt.setp(axs[0].get_xticklabels(), rotation=45)
plt.setp(axs[0].get_yticklabels(), rotation=45)
cbarax1 = fig.add_axes([0.48, 0.2, 0.01, 0.6])
fig.colorbar(demobj, cax=cbarax1, orientation='vertical')
levels = np.linspace(np.min(dataproc[:,-1]),41.93*2-np.min(dataproc[:,-1]),101)
cbobj2 = axs[1].contourf(gravstn_x, gravstn_y, dataproc[:,-1].reshape(np.flip(stn_array_size)),levels=levels,cmap='bwr')
for c in cbobj2.collections:
c.set_edgecolor("face")
axs[1].set_title(r'$\beta$ ($\mu$Gal/mH$_2$O)')
axs[1].set_aspect('equal')
axs[1].grid(c='k')
plt.gca().invert_yaxis()
plt.setp(axs[1].get_xticklabels(), rotation=45)
cbarax2 = fig.add_axes([0.91, 0.2, 0.01, 0.6])
fig.colorbar(cbobj2, cax=cbarax2, orientation='vertical')
return dataproc