Expose CoCiP ellipse inclination

[aaronsarna]

CoCiP models the contrail as an ellipse that may be rotated with respect to the vertical and horizontal axes. The width and depth parameters that pycontrails outputs appear to be the projections of this ellipse onto the vertical and horizontal axes (like B and D in Fig 1 of the CoCiP paper), and as far as I can tell, the angle of inclination is not directly exposed anywhere in pycontrails outputs.

I thought maybe it could be backed out using a combination of width, depth, sigma_yz, and maybe area_eff, cos_a, and sin_a, but conversations with @zebengberg and @roger-teoh indicated that perhaps that isn't sufficient, or at least you might end up missing a sign of the angle of inclination.

Having this information is important when reprojecting CoCiP outputs for alignment with observations in off-nadir geostationary imagery and ground cameras. Could this information be exposed?

(also tagging @thabbott)

[thabbott]

Not sure if this is a solution you talked about with @zebengberg and @roger-teoh, but I think you should be able to get inclination angles by computing eigenvectors of the plume covariance matrix:

import matplotlib.pyplot as plt
import numpy as np

width = 2.0
depth = 1.0
sigma_yz = 0.2

sigma_yy = width**2/8  # schumann 2012 eq. 8
sigma_zz = depth**2/8  # schumann 2012 eq. 8

cov = np.array([sigma_yy, sigma_yz, sigma_yz, sigma_zz]).reshape((2,2))
eigval, eigvec = np.linalg.eig(cov)
imajor = np.argmax(eigval)
iminor = np.argmin(eigval)

y = np.linspace(-3, 3)
z = np.linspace(-3, 3)
yy, zz = np.meshgrid(y, z)

A = 2*np.pi*np.sqrt(np.linalg.det(cov))  # schumann 2012 eq. 3
c = 1/A*np.exp(-1/(2*np.linalg.det(cov))*(yy**2*sigma_zz - 2*yy*zz*sigma_yz + zz**2*sigma_yy))  # schumann 2012 eq. 1

plt.figure()
plt.contourf(y, z, c, levels=100, cmap="Blues_r")
plt.colorbar()
plt.xlabel("y")
plt.ylabel("z")
lam, v = eigval[imajor], eigvec[:,imajor]
plt.plot([0, lam*v[0]], [0, lam*v[1]], "k-", label="semimajor axis")
lam, v = eigval[iminor], eigvec[:,iminor]
plt.plot([0, lam*v[0]], [0, lam*v[1]], "r-", label="semiminor axis")
plt.legend(loc="upper left")
plt.show()

(There's a little trigonometry I haven't done to get inclination from the eigenvectors, but that should be straightforward.)

I'll talk to Roger and Zeb about exposing something about inclination in pycontrails, but thought I'd share this in case it's useful in the meantime.

[aaronsarna]

Thanks, yeah I think this helps.

This makes me realize that tau_contrail also assumes (I think) that the photons are traveling along the z axis. If I'm considering detectability from a geostationary satellite or a ground camera, I guess I need to recompute the width and depth using the ray from the camera to the contrail, rather than the z-axis. It looks like if I did that I could recompute tau using contrail_properties.contrail_optical_depth with the new width. Does that seem right? It's not clear to me why contrail_optical_depth doesn't take depth into account, unless it's somehow implicit from one of the other inputs.

[thabbott]

Yes, good point--tau_contrail is the optical depth along the z-axis, since that's what the radiative forcing calculations (which involve vertical fluxes) require.

contrail_optical_depth needs ice crystals per unit horizontal area (which will tend to be bigger for a deeper contrail) and computes it from n_ice_per_m (number of ice crystals per unit contrail length) and the contrail width. You could get an equivalent value for optical depth along a ray from a camera by using contrail_optical_depth with the width normal to the new axis, like you suggested.

One thing to be careful about is that the quantity computed by contrail_optical_depth is really just a scale for the contrail optical depth. The actual optical depth varies as you move from the center of the contrail toward the edges. You can use equation 63 in the Schumann paper to get contrail optical depth along the vertical axis as a function of transverse distance. I think you probably want something similar, but in a rotated coordinate system that's aligned with the camera ray. (Or maybe you've worked through this part of the problem already, and really just need a value for tau_contrail along the camera ray, in which case the approach you suggested sounds right to me.)

[aaronsarna]

Yep, I'm using appendix A12 of Schumann 2012 for rasterization, which is essentially just an implementation of equation 63.

Does "ice crystals per unit horizontal area" need to be recomputed if the horizontal plane (which I assume really means the plane normal to the vertical axis) is changing to be relative to a satellite ray?

[thabbott]

Ok, great. Yes, if your rasterization scheme needs tau_contrail along a camera ray then you need ice crystals per unit area normal to the camera ray, which you can get by setting the width value passed to contrail_optical_depth to the width relative to the camera ray.

[aaronsarna]

But the n_ice_per_m would not change if the plane changes?

[thabbott]

Ah, I see. No, it does not change. (It's the number of ice crystals per unit contrail length, so it doesn't depend on the angle you're viewing the contrail from.)

[aaronsarna]

Ok I think I have what I need then. Should I close this issue, or is it worth leaving this open to see if pycontrails wants to implement some of this computation itself?

Also, another thing I had to sort out here for the purposes of this reprojection are the x (heading) and y (horizontal, normal to x) axes of the contrail. I'm generally doing all of the reprojection in cartesian coordinates, and I assume there should be some way to back these out from the sin_a and cos_a outputs, but I had some trouble figuring that out. What I ended up doing is locating each waypoint in cartesian space, subtracting waypoint[n]'s coords from waypoint[n+1]'s coords and orthogonalizing the result wrt z to get x, and then y is just np.cross(x, z). If there's a better way to do this using sin_a and cos_a that could be nice. Otherwise, maybe this is another thing that pycontrails could do for you.

[thabbott]

No, let's leave this issue open for now. I might add some orientation outputs when I have spare time, but if you have what you need I won't make it a super high priority.

What you're doing sounds reasonable to me. I'm not sure exactly how you're defining your cartesian coordinate system(s), but if you want to work on a tangent plane at the location of each contrail segment then the unit vector along the x axis of the contrail is just cos_a \hat{x} + sin_a \hat{y}, where \hat{x} and \hat{y} are unit vectors in the eastward and northward directions. My guess is you already know that, though, and for some reason that's not very useful?

[aaronsarna]

I guess the trickiness is whether the tangent plane is tangent to the contrail segment, which may not be parallel to the earth's surface, or instead to the ~sphere of the earth, inflated to the contrail segment's altitude. For a short enough segment it probably doesn't make a difference, though. The documentation could generally be improved here, since I don't think it discusses a tangent plane at all.

FWIW, the 3D cartesian coordinate system I'm using has the origin at the center of the earth, and positive axes pointing to lat/lngs (0, 0), (0, 90) and the north pole. Defining the tangent plane this way is pretty straightforward, but determining what "north" and "east" mean in that context feels a bit less straightforward. Certainly doable, but the approach I ended up taking seemed simpler, at least to me.

[thabbott]

By "tangent plane" I meant a plane tangent to the surface. (If you really want to be careful I think it's tangent to the surface at the location of the earlier of the two waypoints that define a contrail segment.) Agreed that the documentation for contrail outputs is not as detailed as it could be....

It makes sense that you'd want to use a single (global) cartesian coordinate system for your rasterization, though. I think in that case what you're doing (getting segment orientation from waypoint coordinates in your coordinate system) seems simpler than trying to go from spherical coordinates to local cartesian coordinates and then to global cartesian coordinates.