/
sphere_location.py
160 lines (135 loc) · 6.93 KB
/
sphere_location.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
import mcm_utils
import numpy as np
import math
import scipy
from mcm_utils import deg2rad, rad2deg
from scipy.interpolate import CloughTocher2DInterpolator
from scipy.interpolate import NearestNDInterpolator
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import phys_utils
from mpl_toolkits.basemap import Basemap
class sphere_coordinates:
def __init__(self,segments=30):
self.segments = segments
self.build_sphere_mapping()
def build_sphere_mapping(self):
longitudes=np.linspace(0,360,num=self.segments,endpoint=False)
extended_longitudes=np.tile(longitudes,(self.segments,1))
sin_latitudes=np.linspace(1,-1,num=self.segments,endpoint=False)
self.longitudes=np.linspace(-180,180,num=self.segments,endpoint=True)
self.latitudes=rad2deg(np.arcsin(np.linspace(1,-1,num=self.segments,endpoint=True)))
extended_sin_latitudes=np.tile(sin_latitudes,(self.segments,1)).T
self.upper_right=np.swapaxes(np.array([extended_longitudes,extended_sin_latitudes]),0,2)
self.rights=extended_longitudes
self.uppers=extended_sin_latitudes
def direction_to_region(self,direction_vectors):
latlongs = mcm_utils.geographic_coordinates(direction_vectors).T
return self.lat_longs_to_region(latlongs)
def lat_longs_to_region(self,lat_longs):
xreg=((lat_longs[1]+180)/360*self.segments).astype(int)
yreg=((1-np.sin(mcm_utils.deg2rad(lat_longs[0])))*(self.segments/2)).astype(int)
return [yreg, xreg]
class geocentric_data:
def convert_lambert(self):
lambert_coordinate=np.sin(deg2rad(self.latitude))
return np.array([self.longitude, lambert_coordinate])
#data must be a sorted array of latitudes, longitudes, and then grid mesh of the coordinates.
def __init__(self,latitudes, longitudes, data, radius=200, discrete=False):
self.radius=radius
self.shape=data.shape
self.latitude=latitudes
self.longitude=longitudes
self.values=data
self.discrete=discrete
self.lambert_data=self.convert_lambert()
self.lambert_x=self.lambert_data[0]
self.lambert_y=self.lambert_data[1]
coordinates=np.swapaxes(np.array(np.meshgrid(self.latitude,self.longitude)),0,2)
interpshape=(coordinates.shape[0]**2,2)
flatshape=(coordinates.shape[0]**2,)
cartshape=(coordinates.shape[0]**2,3)
if self.discrete:
self.function_interpolator=CloughTocher2DInterpolator(np.reshape(coordinates,interpshape),np.reshape(self.values,flatshape))
else:
self.function_interpolator=NearestNDInterpolator(np.reshape(coordinates,interpshape),np.reshape(self.values,flatshape))
if not discrete:
r=self.values
phi=deg2rad(self.longitude)
theta=deg2rad(90-self.latitude)
sphere_coordinates=np.meshgrid(theta,phi)
theta_mesh=sphere_coordinates[0]
phi_mesh=sphere_coordinates[1]
sphere_coordinates=np.swapaxes(sphere_coordinates,0,2)
sphere_gradient=np.gradient(self.values,math.pi/self.latitude.shape[0],(2*math.pi)/self.longitude.shape[0])
drdtheta=np.squeeze(sphere_gradient[0])
drdphi=np.squeeze(sphere_gradient[1])
cos_theta=np.cos(theta_mesh)
sin_theta=np.sin(theta_mesh)
cos_phi=np.cos(phi_mesh)
sin_phi=np.sin(phi_mesh)
dzdr=cos_theta
dzdtheta=drdtheta*cos_theta-r*sin_theta
dzdphi=drdphi*cos_theta
dzdx=cos_theta*sin_phi*dzdr-(sin_theta/(r*sin_phi))*dzdtheta+(cos_theta*cos_phi/r)*dzdphi
dzdy=sin_theta*sin_phi*dzdr+(cos_theta/(r*sin_phi))*dzdtheta+(sin_theta*cos_phi/r)*dzdphi
gradients=np.array([dzdx,dzdy])
self.gradient_interpolator=CloughTocher2DInterpolator(np.reshape(coordinates,interpshape),np.reshape(gradients,interpshape))
def visualize_lambert(self,mapview=False,log_scale=True,show=True,*args,**vargs):
if mapview:
# llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon
# are the lat/lon values of the lower left and upper right corners
# of the map.
# resolution = 'c' means use crude resolution coastlines.
m = Basemap(projection='cea',llcrnrlat=-90,urcrnrlat=90,llcrnrlon=-180,urcrnrlon=180,resolution='c')
mapx,mapy=m(self.longitude,self.latitude)
#print(self.values)
m.drawcoastlines(linewidth=0.25)
m.pcolormesh(mapx,mapy,self.values,
norm=colors.SymLogNorm(linthresh=0.03, linscale=0.03,
vmin=self.values.min(), vmax=self.values.max()),*args,**vargs)
#m.fillcontinents(color='coral',lake_color='aqua')
# print(mapx)
# draw parallels and meridians.
#m.drawparallels(np.arange(-90.,91.,30.))
#m.drawmeridians(np.arange(-180.,181.,60.))
#m.drawmapboundary(fill_color='aqua')
else:
plt.pcolormesh(self.lambert_x,self.lambert_y,self.values,cmap=cmap,
norm=colors.SymLogNorm(linthresh=0.03, linscale=0.03,
vmin=self.values.min(), vmax=self.values.max()),*args,**vargs)
if show:
plt.show()
return plt
def interpolate_data(self, vectors, system='cartesian'):
if system=='cartesian':
return self.function_interpolator(mcm_utils.geographic_coordinates(vectors))
if system=='spherical':
latlongs=np.array([90-rad2deg(vectors[2]),rad2deg(vectors[1])]).T
return self.function_interpolator(latlongs)
if system=='geographic':
return self.function_interpolator(vectors)
def interpolate_gradient(self,vectors, system='cartesian'):
if self.discrete:
return None
if system=='cartesian':
return self.gradient_interpolator(mcm_utils.geographic_coordinates(vectors))
if system=='spherical':
latlongs=np.array([90-rad2deg(vectors[2]),rad2deg(vectors[1])]).T
return self.gradient_interpolator(latlongs)
if system=='geographic':
return self.gradient_interpolator(vectors)
if __name__=='__main__':
longitudes=np.linspace(-180,180,num=20,endpoint=False)
latitudes=np.linspace(-90,90,num=20,endpoint=False)
#data=np.zeros((20,20))
#for i in range(latitudes.shape[0]):
# for k in range(longitudes.shape[0]):
# #print("for latitude",latitudes[i],"and for longitude",longitudes[k])
# profile=phys_utils.electron_density_profile(latitudes[i],longitudes[k])
# data[i][k]=np.amax(profile)
#np.save('electron_density',data)
data=np.load('electron_density.npy')
constant=geocentric_data(latitudes,longitudes,data)
maxvector=np.array([20,50])
# print(constant.interpolate_gradient(maxvector,system='geographic'))