"""Planetary mapping routines.
phi = 0 faces toward the observer
phi = pi thus faces away from the observer
theta=pi/2 is the z-axis or 'north pole'
theta=-pi/2 is the 'south pole'
"""
# 2010-01-15 20:31 IJC: Started . . .
# 2013-08-07 11:04 IJMC: Added mu field for cells and maps
from numpy import pi
import numpy as np
[docs]def polyarea(x, y):
"""Compute the area of a polygon whose vertices are at the points (x,y).
:INPUTS:
x, y : 1D sequences
Cartesian coordinates of the (non-intersecting) polygon.
:REFERENCE:
http://mathworld.wolfram.com/PolygonArea.html
"""
# 2013-05-29 12:18 IJMC: Created
area = 0.
npts = max(len(x), len(y))
for ii in xrange(npts):
area += x[ii]*y[(ii+1) % npts] - x[(ii+1) % npts]*y[ii]
return np.abs(area*0.5)
[docs]def makegrid(nphi,ntheta):
"""Make grids of phi and theta values with the specified number of
points in each direction. Phi ranges from 0 to 2pi, and theta
ranges from -pi/2 to pi/2.
Returns (phi, theta)"""
# 2010-01-15 20:29 IJC: Created
# 2013-08-18 15:57 IJMC: Updated so phi values don't repeat at 0 & 2pi
from numpy import meshgrid, linspace
phi,theta = meshgrid(linspace(0,2*pi,nphi+1)[0:-1],linspace(-pi/2,pi/2,ntheta))
return phi,theta
[docs]def rotcoord(phi,theta,iangle,rangle):
"""rotate coordinate system from local (planetary) coordinates into
observer-oriented coordinages. rangle and iangle are the rotation
and inclination angles of the planet in radians.
returns (phi2,theta2)
phi2 will be in the range (0,2pi) and
theta2 will be in the range (-pi/2,pi/2)
"""
# 2010-01-15 21:09 IJC: Created
from numpy import array,cos,sin,vstack,dot,arctan2,sqrt,pi, abs
nphi,ntheta = phi.shape
phir = phi+rangle
costheta = cos(theta)
x = cos(phir)*costheta
y = sin(phir)*costheta
z = sin(theta)
xyz = vstack((x.ravel(),y.ravel(),z.ravel()))
beta = pi/2.-iangle
sinbeta = sin(beta)
cosbeta = cos(beta)
rotmat = array([[cosbeta,0,-sinbeta],[0,1,0],[sinbeta,0,cosbeta]])
xyz2 = dot(rotmat,xyz)
x2,y2,z2 = xyz2[0].reshape(nphi,ntheta),xyz2[1].reshape(nphi,ntheta),xyz2[2].reshape(nphi,ntheta)
phi2 = arctan2(y2,x2) % (2*pi)
theta2 = arctan2(z2,sqrt(x2**2+y2**2))
return phi2,theta2
[docs]def visiblemap(phi,theta,iangle,rangle):
"""Return a 2D boolean map that's True for the planetary latitude,
longitude values visible from an observer. For rangle (rotation
angle) zero and iangle (inclination angle) equal to pi/2, phi=pi
is toward the observer.
"""
# 2010-01-18 08:13 IJC: Created
nphi,ntheta = phi.shape
phi2,theta2 = rotcoord(phi,theta,iangle,rangle)
vismap = (phi2<pi/2)+(phi2>1.5*pi)
return vismap
[docs]def wedgemap(phi,theta,phi0,phi1):
"""Return a 2D boolean map that's True for the latitude,
longitude values in a given longitudinal planet 'wedge'.
phi0,phi1 should be in the range (0,2*pi)
"""
# 2010-01-18 08:13 IJC: Created
phi0 = phi0 % (2*pi)
phi1 = phi1 % (2*pi)
if phi1>=phi0:
wedge = (phi<phi1)*(phi>=phi0)
elif phi1<phi0:
wedge = (phi>=phi0)+(phi<phi1)
return wedge
[docs]def wedgestack(phi,theta,nwedge,phi0):
"""Return a stack of 2D boolean maps via wedgemap.
phi0 defines the center of wedge zero.
"""
# 2010-01-18 08:13 IJC: Created
# 2013-12-16 05:39 IJMC: Slight speed boost.
from numpy import zeros,arange
dphi = 2*pi/nwedge
pp0=arange(-.5,nwedge-.5,dtype=float)*dphi+phi0
pp1=pp0 + dphi #arange(.5,nwedge+.5,dtype=float)*dphi+phi0
wedges = zeros((nwedge,)+phi.shape,float)
for ii in range(nwedge):
wedges[ii,:,:]= wedgemap(phi,theta,pp0[ii],pp1[ii])
return wedges
[docs]def projarea(phi,theta):
"""Return a 2D map of the projected area. Note that you need to
determine for yourself which parts of the areal map are actually
visible to the observer.
Assumes phi and theta are grids straight out of meshgrid.
Assumes dA = cos(t) dt df (theta=t, phi=f)
and thus da = cos(f) cos(t)**2 dt df
EXAMPLE:
import maps, pylab
phi, theta = maps.makegrid(300,240)
vis = maps.visiblemap(phi,theta,pi/2,0)
da = maps.projarea(phi,theta)
pylab.imshow(da*vis)
pylab.title('visible projected area sums to: %f (pi)' % (da*vis).sum())
"""
# 2010-01-18 09:23 IJC: Created
# 2013-08-18 21:23 IJMC: Updated documentation, used median instead.
from numpy import cos
df = np.median(np.diff(phi[0])) #phi[0,1]-phi[0,0]
dt = np.median(np.diff(theta[:,0])) #theta[1,0]-theta[0,0]
da = (cos(theta)**2) * cos(phi) * dt * df
return da
[docs]def wedgebasis(phi,theta,iangle,rangle,nwedge,phi0, fwedge=None):
"""Return a set of basis functions for the flux from each wedge.
phi,theta are the observer-centered grids from MAPS.MAKEGRID
iangle is the inclination of the system (0 = pole-on)
rangle (seq.) is a sequence of (0,2pi) rotation values at which
the wedge-based flux is evaluated.
nwedge is the number of wedges
phi0 defines the center of wedge zero.
fwedge -- an optional sequence to set individual flux values for each wedge.
"""
# 2010-01-18 08:13 IJC: Created
from numpy import array,arange,zeros,tile
rangle = array(rangle,copy=True)
nrot = len(rangle)
bases = zeros((nrot,nwedge),float)
da = projarea(phi,theta)
vis = (phi<=0.5*pi) + (phi>1.5*pi)
davis = da*vis
phi2,theta2 = rotcoord(phi,theta,iangle,rangle[0])
for ii in range(nrot):
# Faster to shift reference point than entire phi2 grid:
wedges = wedgestack(phi2,theta2,nwedge,phi0-rangle[ii]+rangle[0])
bases[ii,:] = (wedges*davis).sum(2).sum(1)
if fwedge<>None:
fwedge = tile(fwedge,(nrot,1))
bases = bases * fwedge
return bases
[docs]def errseries(param,phi,theta,rangle,meas,err, retmodel=False):
"""Give the chi-squared to an input timeseries to get the wedge
coefficients and inclination angle.
param -- [inclination, phi0, wedge coefficients]
phi,theta -- from MAKEGRID
rangle -- rotation angles at which flux was measured
meas -- measured flux
err -- error on measurement (one-sigma)
retmodel -- whether to return the model, instead of the chi-squared.
Usage:
from scipy import optimize
optimize.fmin(maps.errseries, guess, args=(...))
"""
# 2010-01-18 14:04 IJC: Created
# 2013-12-16 05:17 IJMC: Added 'retmodel' option.
iangle = param[0]
phi0 = param[1]
fcoef = param[2::]
nwedge = len(fcoef)
model = wedgebasis(phi,theta,iangle,rangle,nwedge,phi0,fwedge=fcoef).sum(1)
if retmodel:
ret = model
else:
chisq = (((model-meas)/err)**2).sum()
print "param, chisq>>" +str(param)+', '+str(chisq)
ret = chisq
return ret
[docs]def makespot(spotlat, spotlon, spotrad, phi, theta):
"""
:INPUTS:
spotlat : scalar
Latitude of spot center, in radians, from 0 to pi
spotlon : scalar
Longitude of spot center, in radians, from 0 to 2pi
spotrad : scalar
Radius of spot, in radians.
phi, theta : 2D NumPy arrays
output from :func:`makegrid`. Theta ranges from -pi/2 to +pi/2.
:EXAMPLE:
::
import maps
nlat, nlon = 60, 30
phi, theta = maps.makegrid(nlat, nlon)
# Make a small spot centered near, but not at, the equator:
equator_spot = maps.makespot(0, 0, 0.4, phi, theta)
# Make a larger spot centered near, but not at, the pole:
pole_spot = maps.makespot(1.2, 0, 0.7, phi, theta)
::
import maps
nlat, nlon = 60, 30
map = maps.map(nlat, nlon, i=0., deltaphi=0.)
phi = map.corners_latlon.mean(2)[:,1].reshape(nlon, nlat)
theta = map.corners_latlon.mean(2)[:,0].reshape(nlon, nlat) - np.pi/2.
# Make a small spot centered near, but not at, the equator:
equator_spot = maps.makespot(0, 0, 0.4, phi, theta)
# Make a larger spot centered near, but not at, the pole:
pole_spot = maps.makespot(1.2, 0, 0.7, phi, theta)
"""
# 2013-08-18 16:01 IJMC: Created
pi2 = 0.5*np.pi
xyz = np.array((np.cos(phi) * np.sin(theta + pi2), np.sin(phi) * np.sin(theta + pi2), np.cos(theta + pi2))).reshape(3, phi.size)
# First rotate around z axis, to align spot with sub-observer meridian
# Then, rotate around y axis, to align spot with pole.
zrot = np.array([[np.cos(np.pi-spotlon), -np.sin(np.pi-spotlon), 0], [np.sin(np.pi-spotlon), np.cos(np.pi-spotlon), 0.], [0,0,1]])
yrot = np.array([[np.cos(spotlat+pi2), 0, np.sin(spotlat+pi2)], [0,1,0], [-np.sin(spotlat+pi2), 0, np.cos(spotlat+pi2)]])
xyz = np.dot(np.dot(yrot, zrot), xyz)
# Convert Cartesian to spherical coordinates
ang = np.arccos(xyz[2])
# Spot is where (theta - theta_pole) < radius.
spotmap = ang.T <= spotrad
return spotmap.reshape(phi.shape)
[docs]def makespot_old(phi, theta, da=None, vmap=None, inc=None, rot=None, long=0, lat=0, siz=pi/4, plotalot=False):
"""Make a spot-map.
EXAMPLE:
import maps
inc, rot = pi/2., 0.
phi, theta = maps.makegrid(120,60)
phi2, theta2 = maps.rotcoord(phi,theta,inc,rot)
vmap = maps.visiblemap(phi,theta,inc,rot)
maps.makespot(phi,theta,inc=pi/2,rot=0.,lat=pi/4,long=1.*pi)
maps.makespot(phi2,theta2,vmap=vmap,lat=pi/4,long=1.*pi)
"""
if plotalot:
from pylab import *
from tools import nextfig
figure(nextfig(), [15,15]); clf()
if da is None:
da = projarea(phi,theta)
if plotalot:
subplot(3,3,1)
imshow(phi); colorbar(); title('phi')
subplot(3,3,2)
imshow(theta); colorbar(); title('theta')
subplot(3,3,3)
imshow(da); colorbar(); title('d_area')
if inc is not None and rot is not None:
phi, theta = rotcoord(phi,theta,inc,rot)
if vmap is None:
vmap = visiblemap(phi,theta,inc,rot)
if plotalot:
subplot(3,3,4)
imshow(phi); colorbar(); title('phi')
subplot(3,3,5)
imshow(theta); colorbar(); title('theta')
subplot(3,3,6)
imshow(vmap); title('vmap: inc=%s, rot=%s'%(inc,rot)); colorbar()
if long is not None and lat is not None:
phi, theta, = rotcoord(phi,theta,-lat, -long)
spot = vmap*da*((theta-theta.max())>-siz)
if plotalot:
subplot(3,3,7)
imshow(phi); colorbar(); title('phi')
subplot(3,3,8)
imshow(theta); colorbar(); title('theta')
subplot(3,3,9)
imshow(spot); title('spot: lat=%s, long=%s'%(lat,long)); colorbar()
return phi, theta, da, vmap, spot, (inc, rot, long, lat, siz)
class mapcell:
def __init__(self):
self.corners = np.zeros((3, 4), dtype=float)
self.corners_latlon = np.zeros((2, 4), dtype=float)
self.vcorners = np.zeros((3, 4), dtype=float)
self.rvcorners = np.zeros(4, dtype=float)
self.visible_corners = np.zeros((3, 4), dtype=float)
self.visible_vcorners = np.zeros((3, 4), dtype=float)
self.visible_rvcorners = np.zeros(4, dtype=float)
self.projected_area = 0.
self.mu = 0.
return
def get_mu(self):
### Compute mu:
normal_vector = np.dot(np.linalg.pinv(self.corners.T), np.ones(4))
self.mu = normal_vector[0] / np.sqrt(np.dot(normal_vector, normal_vector))
return
def get_projected_area(self, i):
if (self.corners[0] <= 0).all():
# cell is hidden, on the back side.
area = 0.
self.visible_corners = self.corners * np.nan
elif (self.corners[0] > 0).all():
# cell is completely visible, on the front side.
self.visible_corners = self.corners
y = self.corners[1]
z = self.corners[2]
inds = np.argsort(np.arctan2(z-z.mean(), y-y.mean()))
area = polyarea(y[inds], z[inds])
else:
# Cell is only partially visible (on the limb). Find the
# nearest point on on the limb, with the same latitude as
# each vertex.
visible_corners = self.corners.copy()
back_indices = (visible_corners[0] < 0).nonzero()[0]
for ii in back_indices:
newx = 0. # on the limb!
newy = np.sin(self.corners_latlon[0,ii]) * \
np.sqrt(1. - np.tan(i)**2 / np.tan(self.corners_latlon[0,ii])**2)
if visible_corners[1,ii]/newy < 0:
newy *= -1
newz = np.cos(self.corners_latlon[0,ii]) / np.cos(i)
visible_corners[:, ii] = newx, newy, newz
if not (np.isfinite(visible_corners)).all():
self.visible_corners = self.corners * np.nan
area = 0
print "Non-finite projected corners; need to fix this."
else:
self.visible_corners = visible_corners
y = self.visible_corners[1]
z = self.visible_corners[2]
yz = np.unique(zip(y,z))
inds = np.argsort(np.arctan2(yz[:,1]-yz[:,1].mean(), yz[:,0]-yz[:,0].mean()))
area = polyarea(yz[inds,0], yz[inds,1])
self.projected_area = area
return
[docs]class map:
"""Very handy spherical mapping object.
:INPUTS:
nlon, nlat : scalars
If mod=='latlon', these inputs specify the number of grid cells
across map, in latitude and longitude.
i : scalar
the inclination, is in units of radians. Zero means we see the
object equator-on; pi/2 means we see it pole-on.
type : str
For now, only valid entry is 'latlon'. Eventually, this could be
expanded to allow other types of projection grids.
deltaphi : scalar
Rotation of map, specified in radians.
:OUTPUT:
A map-class object with various useful fields. Most of these
fields refer to the coordinates (either Cartesian or spherical
polar) or the projected radial velocities at the corners of
specified grid cells, or the approximate projected areas of
these grid cells.
:NOTES:
I have *not* been as careful as I should be in this code -- my
original goal was speed rather than exactitude. This means that
some values are returned as 'nan', and the projected areas are
only roughly correct. There's plenty of room for improvement!
"""
# 2013-05-29 09:37 IJMC: Created
# 2013-08-07 11:05 IJMC: Added mu field for maps and cells
# 2014-08-07 15:00 IJMC: Updated documentation -- exactly 1 year later!
def __init__(self, nlon=20, nlat=10, type='latlon', deltaphi=0, i=0):
self.type = type
self.nlon = nlon
self.nlat = nlat
self.ncell = nlon*nlat
self.deltaphi = deltaphi
self.i = i
self.cells = []
self.visible_corners = np.zeros((self.ncell, 3, 4), dtype=float)
self.corners = np.zeros((self.ncell, 3, 4), dtype=float)
self.corners_latlon = np.zeros((self.ncell, 2, 4), dtype=float)
self.rvcorners = np.zeros((self.ncell, 4), dtype=float)
self.visible_rvcorners = np.zeros((self.ncell, 4), dtype=float)
self.projected_area = np.zeros(self.ncell, dtype=float)
self.mu = np.zeros(self.ncell, dtype=float)
self.phi = np.zeros(self.ncell)
self.theta = np.zeros(self.ncell)
### Initialize coordinate system:
phi0 = np.arange(0, self.nlon+1) * (2*np.pi/self.nlon)
theta0 = np.arange(0, self.nlat+1) * (np.pi/self.nlat)
phi, theta = np.meshgrid(phi0, theta0)
### Rotate by deltaPhi:
phi1 = (phi + deltaphi).ravel()
theta1 = theta.ravel()
### Convert to x1, y1, z1:
xyz1 = np.vstack((np.sin(theta1) * np.cos(phi1), \
np.sin(theta1) * np.sin(phi1), \
np.cos(theta1)))
### Rotate by inclination angle i:
rot_matrix = np.array([[np.cos(i), 0, -np.sin(i)], \
[0,1,0], [np.sin(i), 0, np.cos(i)]])
xyz2 = np.dot(rot_matrix, xyz1)
xyz3 = xyz2.reshape(3, nlat+1, nlon+1)
kk = 0
for ii in xrange(self.nlat):
for jj in xrange(self.nlon):
cell = mapcell()
cell.corners = xyz3[:, ii:ii+2, jj:jj+2].reshape(3,4)
cell.corners_latlon = np.vstack((theta[ii:ii+2,jj:jj+2].ravel(), phi[ii:ii+2,jj:jj+2].ravel()))
cell.rvcorners = xyz3[1,ii:ii+2,jj:jj+2].ravel() * np.cos(i)
cell.get_projected_area(i)
cell.get_mu()
cell.visible_rvcorners = cell.visible_corners[1] * np.cos(i)
self.cells.append(cell)
self.corners[kk] = cell.corners
self.visible_corners[kk] = cell.visible_corners
self.projected_area[kk] = cell.projected_area
self.mu[kk] = cell.mu
self.corners_latlon[kk] = cell.corners_latlon
self.rvcorners[kk] = cell.rvcorners
self.visible_rvcorners[kk] = cell.visible_rvcorners
kk += 1
return None
[docs] def get_vprofile(self, v):
"""Compute velocity profile for normalized velocity values v.
:INPUTS:
v : NumPy array
Velocity normalized by the maximum rotation velocity
observed; i.e., to convert v to true velocities, multiply
by 2piR/P.
"""
# 2013-05-29 12:28 IJMC: Created
profile = np.zeros(v.shape, dtype=float)
for ii in xrange(self.ncell):
vmin, vmax = self.visible_rvcorners[ii].min(), self.visible_rvcorners[ii].max()
profile[(v > vmin) * (v <= vmax)] += 1 #self.projected_area[ii]
return profile