def plotArc(start_angle, stop_angle, radius, width, **kwargs):
""" write a docstring for this function"""
numsegments = 100
theta = np.radians(np.linspace(start_angle+90, stop_angle+90, numsegments))
centerx = 0
centery = 0
x1 = -np.cos(theta) * (radius)
y1 = np.sin(theta) * (radius)
stack1 = np.column_stack([x1, y1])
x2 = -np.cos(theta) * (radius + width)
y2 = np.sin(theta) * (radius + width)
stack2 = np.column_stack([np.flip(x2, axis=0), np.flip(y2,axis=0)])
#add the first values from the first set to close the polygon
np.append(stack2, [[x1[0],y1[0]]], axis=0)
arcArray = np.concatenate((stack1,stack2), axis=0)
return patches.Polygon(arcArray, True, **kwargs), ((x1, y1), (x2, y2))
python类sin()的实例源码
def ct2lg(dX, dY, dZ, lat, lon):
n = dX.size
R = np.zeros((3, 3, n))
R[0, 0, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.cos(np.deg2rad(lon)))
R[0, 1, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.sin(np.deg2rad(lon)))
R[0, 2, :] = np.cos(np.deg2rad(lat))
R[1, 0, :] = -np.sin(np.deg2rad(lon))
R[1, 1, :] = np.cos(np.deg2rad(lon))
R[1, 2, :] = np.zeros((1, n))
R[2, 0, :] = np.multiply(np.cos(np.deg2rad(lat)), np.cos(np.deg2rad(lon)))
R[2, 1, :] = np.multiply(np.cos(np.deg2rad(lat)), np.sin(np.deg2rad(lon)))
R[2, 2, :] = np.sin(np.deg2rad(lat))
dxdydz = np.column_stack((np.column_stack((dX, dY)), dZ))
RR = np.reshape(R[0, :, :], (3, n))
dx = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
RR = np.reshape(R[1, :, :], (3, n))
dy = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
RR = np.reshape(R[2, :, :], (3, n))
dz = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
return dx, dy, dz
def distance(self, lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6371 * c
return km
def distance(self, lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = numpy.sin(dlat/2)**2 + numpy.cos(lat1) * numpy.cos(lat2) * numpy.sin(dlon/2)**2
c = 2 * numpy.arcsin(numpy.sqrt(a))
km = 6371 * c
return km
def ct2lg(dX, dY, dZ, lat, lon):
n = dX.size
R = np.zeros((3, 3, n))
R[0, 0, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.cos(np.deg2rad(lon)))
R[0, 1, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.sin(np.deg2rad(lon)))
R[0, 2, :] = np.cos(np.deg2rad(lat))
R[1, 0, :] = -np.sin(np.deg2rad(lon))
R[1, 1, :] = np.cos(np.deg2rad(lon))
R[1, 2, :] = np.zeros((1, n))
R[2, 0, :] = np.multiply(np.cos(np.deg2rad(lat)), np.cos(np.deg2rad(lon)))
R[2, 1, :] = np.multiply(np.cos(np.deg2rad(lat)), np.sin(np.deg2rad(lon)))
R[2, 2, :] = np.sin(np.deg2rad(lat))
dxdydz = np.column_stack((np.column_stack((dX, dY)), dZ))
RR = np.reshape(R[0, :, :], (3, n))
dx = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
RR = np.reshape(R[1, :, :], (3, n))
dy = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
RR = np.reshape(R[2, :, :], (3, n))
dz = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
return dx, dy, dz
def distance(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6371 * c
return km
def todictionary(self, time_series=False):
# convert the ETM adjustment into a dirtionary
# optionally, output the whole time series as well
# start with the parameters
etm = dict()
etm['Linear'] = {'tref': self.Linear.tref, 'params': self.Linear.values.tolist()}
etm['Jumps'] = [{'type':jump.type, 'year': jump.year, 'a': jump.a.tolist(), 'b': jump.b.tolist(), 'T': jump.T} for jump in self.Jumps.table]
etm['Periodic'] = {'frequencies': self.Periodic.frequencies, 'sin': self.Periodic.sin.tolist(), 'cos': self.Periodic.cos.tolist()}
if time_series:
ts = dict()
ts['t'] = self.ppp_soln.t.tolist()
ts['x'] = self.ppp_soln.x.tolist()
ts['y'] = self.ppp_soln.y.tolist()
ts['z'] = self.ppp_soln.z.tolist()
etm['time_series'] = ts
return etm
def ct2lg(self, dX, dY, dZ, lat, lon):
n = dX.size
R = numpy.zeros((3, 3, n))
R[0, 0, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon)))
R[0, 1, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon)))
R[0, 2, :] = numpy.cos(numpy.deg2rad(lat))
R[1, 0, :] = -numpy.sin(numpy.deg2rad(lon))
R[1, 1, :] = numpy.cos(numpy.deg2rad(lon))
R[1, 2, :] = numpy.zeros((1, n))
R[2, 0, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon)))
R[2, 1, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon)))
R[2, 2, :] = numpy.sin(numpy.deg2rad(lat))
dxdydz = numpy.column_stack((numpy.column_stack((dX, dY)), dZ))
RR = numpy.reshape(R[0, :, :], (3, n))
dx = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
RR = numpy.reshape(R[1, :, :], (3, n))
dy = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
RR = numpy.reshape(R[2, :, :], (3, n))
dz = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
return dx, dy, dz
def distance(self, lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = numpy.sin(dlat/2)**2 + numpy.cos(lat1) * numpy.cos(lat2) * numpy.sin(dlon/2)**2
c = 2 * numpy.arcsin(numpy.sqrt(a))
km = 6371 * c
return km
def distance(self, lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = numpy.sin(dlat/2)**2 + numpy.cos(lat1) * numpy.cos(lat2) * numpy.sin(dlon/2)**2
c = 2 * numpy.arcsin(numpy.sqrt(a))
km = 6371 * c
return km
def ct2lg(self, dX, dY, dZ, lat, lon):
n = dX.size
R = numpy.zeros((3, 3, n))
R[0, 0, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon)))
R[0, 1, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon)))
R[0, 2, :] = numpy.cos(numpy.deg2rad(lat))
R[1, 0, :] = -numpy.sin(numpy.deg2rad(lon))
R[1, 1, :] = numpy.cos(numpy.deg2rad(lon))
R[1, 2, :] = numpy.zeros((1, n))
R[2, 0, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon)))
R[2, 1, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon)))
R[2, 2, :] = numpy.sin(numpy.deg2rad(lat))
dxdydz = numpy.column_stack((numpy.column_stack((dX, dY)), dZ))
RR = numpy.reshape(R[0, :, :], (3, n))
dx = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
RR = numpy.reshape(R[1, :, :], (3, n))
dy = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
RR = numpy.reshape(R[2, :, :], (3, n))
dz = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
return dx, dy, dz
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def plotPlaceTxRxSphereXY(Ax,xtx,ytx,xrx,yrx,x0,y0,a):
Xlim = Ax.get_xlim()
Ylim = Ax.get_ylim()
FS = 20
Ax.scatter(xtx,ytx,s=100,color='k')
Ax.text(xtx-0.75,ytx+1.5,'$\mathbf{Tx}$',fontsize=FS+6)
Ax.scatter(xrx,yrx,s=100,color='k')
Ax.text(xrx-0.75,yrx-4,'$\mathbf{Rx}$',fontsize=FS+6)
xs = x0 + a*np.cos(np.linspace(0,2*np.pi,41))
ys = y0 + a*np.sin(np.linspace(0,2*np.pi,41))
Ax.plot(xs,ys,ls=':',color='k',linewidth=3)
Ax.set_xbound(Xlim)
Ax.set_ybound(Ylim)
return Ax
def calc_IndCurrent_cos_range(self,f,t):
"""Induced current over a range of times"""
Bpx = self.Bpx
Bpz = self.Bpz
a2 = self.a2
azm = np.pi*self.azm/180.
R = self.R
L = self.L
w = 2*np.pi*f
Ax = np.pi*a2**2*np.sin(azm)
Az = np.pi*a2**2*np.cos(azm)
Phi = (Ax*Bpx + Az*Bpz)
phi = np.arctan(R/(w*L))-np.pi # This is the phase and not phase lag
Is = -(w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.cos(w*t + phi)
Ire = -(w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.cos(w*t)*np.cos(phi)
Iim = (w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.sin(w*t)*np.sin(phi)
return Ire,Iim,Is,phi
def is_grid(self, grid, image):
"""
Checks the "gridness" by analyzing the results of a hough transform.
:param grid: binary image
:return: wheter the object in the image might be a grid or not
"""
# - Distance resolution = 1 pixel
# - Angle resolution = 1° degree for high line density
# - Threshold = 144 hough intersections
# 8px digit + 3*2px white + 2*1px border = 16px per cell
# => 144x144 grid
# 144 - minimum number of points on the same line
# (but due to imperfections in the binarized image it's highly
# improbable to detect a 144x144 grid)
lines = cv2.HoughLines(grid, 1, np.pi / 180, 144)
if lines is not None and np.size(lines) >= 20:
lines = lines.reshape((lines.size / 2), 2)
# theta in [0, pi] (theta > pi => rho < 0)
# normalise theta in [-pi, pi] and negatives rho
lines[lines[:, 0] < 0, 1] -= np.pi
lines[lines[:, 0] < 0, 0] *= -1
criteria = (cv2.TERM_CRITERIA_EPS, 0, 0.01)
# split lines into 2 groups to check whether they're perpendicular
if cv2.__version__[0] == '2':
density, clmap, centers = cv2.kmeans(
lines[:, 1], 2, criteria, 5, cv2.KMEANS_RANDOM_CENTERS)
else:
density, clmap, centers = cv2.kmeans(
lines[:, 1], 2, None, criteria,
5, cv2.KMEANS_RANDOM_CENTERS)
if self.debug:
self.save_hough(lines, clmap)
# Overall variance from respective centers
var = density / np.size(clmap)
sin = abs(np.sin(centers[0] - centers[1]))
# It is probably a grid only if:
# - centroids difference is almost a 90° angle (+-15° limit)
# - variance is less than 5° (keeping in mind surface distortions)
return sin > 0.99 and var <= (5*np.pi / 180) ** 2
else:
return False
def save_hough(self, lines, clmap):
"""
:param lines: (rho, theta) pairs
:param clmap: clusters assigned to lines
:return: None
"""
height, width = self.image.shape
ratio = 600. * (self.step+1) / min(height, width)
temp = cv2.resize(self.image, None, fx=ratio, fy=ratio,
interpolation=cv2.INTER_CUBIC)
temp = cv2.cvtColor(temp, cv2.COLOR_GRAY2BGR)
colors = [(0, 127, 255), (255, 0, 127)]
for i in range(0, np.size(lines) / 2):
rho = lines[i, 0]
theta = lines[i, 1]
color = colors[clmap[i, 0]]
if theta < np.pi / 4 or theta > 3 * np.pi / 4:
pt1 = (rho / np.cos(theta), 0)
pt2 = (rho - height * np.sin(theta) / np.cos(theta), height)
else:
pt1 = (0, rho / np.sin(theta))
pt2 = (width, (rho - width * np.cos(theta)) / np.sin(theta))
pt1 = (int(pt1[0]), int(pt1[1]))
pt2 = (int(pt2[0]), int(pt2[1]))
cv2.line(temp, pt1, pt2, color, 5)
self.save2image(temp)
def is_grid(self, grid, image):
"""
Checks the "gridness" by analyzing the results of a hough transform.
:param grid: binary image
:return: wheter the object in the image might be a grid or not
"""
# - Distance resolution = 1 pixel
# - Angle resolution = 1° degree for high line density
# - Threshold = 144 hough intersections
# 8px digit + 3*2px white + 2*1px border = 16px per cell
# => 144x144 grid
# 144 - minimum number of points on the same line
# (but due to imperfections in the binarized image it's highly
# improbable to detect a 144x144 grid)
lines = cv2.HoughLines(grid, 1, np.pi / 180, 144)
if lines is not None and np.size(lines) >= 20:
lines = lines.reshape((lines.size/2), 2)
# theta in [0, pi] (theta > pi => rho < 0)
# normalise theta in [-pi, pi] and negatives rho
lines[lines[:, 0] < 0, 1] -= np.pi
lines[lines[:, 0] < 0, 0] *= -1
criteria = (cv2.TERM_CRITERIA_EPS, 0, 0.01)
# split lines into 2 groups to check whether they're perpendicular
if cv2.__version__[0] == '2':
density, clmap, centers = cv2.kmeans(
lines[:, 1], 2, criteria,
5, cv2.KMEANS_RANDOM_CENTERS)
else:
density, clmap, centers = cv2.kmeans(
lines[:, 1], 2, None, criteria,
5, cv2.KMEANS_RANDOM_CENTERS)
# Overall variance from respective centers
var = density / np.size(clmap)
sin = abs(np.sin(centers[0] - centers[1]))
# It is probably a grid only if:
# - centroids difference is almost a 90° angle (+-15° limit)
# - variance is less than 5° (keeping in mind surface distortions)
return sin > 0.99 and var <= (5*np.pi / 180) ** 2
else:
return False
def build_2D_cov_matrix(sigmax,sigmay,angle,verbose=True):
"""
Build a covariance matrix for a 2D multivariate Gaussian
--- INPUT ---
sigmax Standard deviation of the x-compoent of the multivariate Gaussian
sigmay Standard deviation of the y-compoent of the multivariate Gaussian
angle Angle to rotate matrix by in degrees (clockwise) to populate covariance cross terms
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
covmatrix = tu.build_2D_cov_matrix(3,1,35)
"""
if verbose: print ' - Build 2D covariance matrix with varinaces (x,y)=('+str(sigmax)+','+str(sigmay)+\
') and then rotated '+str(angle)+' degrees'
cov_orig = np.zeros([2,2])
cov_orig[0,0] = sigmay**2.0
cov_orig[1,1] = sigmax**2.0
angle_rad = (180.0-angle) * np.pi/180.0 # The (90-angle) makes sure the same convention as DS9 is used
c, s = np.cos(angle_rad), np.sin(angle_rad)
rotmatrix = np.matrix([[c, -s], [s, c]])
cov_rot = np.dot(np.dot(rotmatrix,cov_orig),np.transpose(rotmatrix)) # performing rot * cov * rot^T
return cov_rot
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def residual(r,theta,u,d):
out = np.empty_like(u)
out[0] = (2*np.sin(theta)*r*d(u[0],1,0)
+ r*r*np.sin(theta)*d(u[0],2,0)
+ np.cos(theta)*d(u[0],0,1)
+ np.sin(theta)*d(u[1],0,2))
out[1] = (2*np.sin(theta)*r*d(u[1],1,0)
+ r*r*np.sin(theta)*d(u[1],2,0)
+ np.cos(theta)*d(u[1],0,1)
+ np.sin(theta)*d(u[1],0,2))
return out
def residual(r,theta,u,d):
u = u[0]
out = (2*np.sin(theta)*r*d(u,1,0)
+ r*r*np.sin(theta)*d(u,2,0)
+ np.cos(theta)*d(u,0,1)
+ np.sin(theta)*d(u,0,2))
out = out.reshape(tuple([1]) + out.shape)
return out
