SHDOM single-scattering adjoint

This folder conatains the most recent (3 May 2015) single-scattering adjoint calculations from frank.

This has the log output files (.log), radiance output files (.arad), and optical property adjoint files (.adj) for the "true" and "guess" cases. The first set of output is for the single scattering forward radiance calculations

les0822nh15t13y135_ocaer1_w0.646_ns1true.*,

les0822nh15t13y135_ocaer1_w0.646_ns1guess.*,

and the second set of outputs is for the full multiple scattering forward radiance calculations,

les0822nh15t13y135_ocaer1_w0.646_ns1true2.*,

les0822nh15t13y135_ocaer1_w0.646_ns1guess2.*.

The cloud field used in these simulations is a 2D slice taken at $y_{(i_y=135)}=8.375 \mathrm{km}$ of the 320x320 high resolution simulation shown below.

The adjoint calculation is for the gradient of the misfit function,

\begin{align} \frac{\partial \Phi(\sigma, \omega)}{\partial(\sigma)} &= \left < \Delta p, \mathcal{U}_{\sigma}[\Delta f] \right > \end{align}

Only the first order of scattering completed so the results are approximate. In the calculation labelled "1" the forward and adjoint calculations are only single-scattering. In the calculation labelled "2", the Forward calculation includes multiple scatterng while the adjoint does not.


In [1]:
%pylab inline

# Imports
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

import pandas as pd


Populating the interactive namespace from numpy and matplotlib

In [ ]:


In [2]:
# Utility functions

def print_head(fname, Nlines=10, indent="\t"):
    "Print the head of the file."
    
    # Print a message and then the first N lines
    print("Showing head: {}".format(fname))
    for i, line in zip(range(Nlines), open(fname, 'r')):
        print(indent + line.strip())
        
    print('\n')

In [3]:
ls


adjoint-derivative-single-scattering.jpg
adjoint-extinction-difference-and residual.jpg
adjoint-extinction-difference-and-residual.jpg
exploring_adjoint.ipynb
exploring_adjoint_with_aerosol.ipynb
les0822nh15t13y135_ocaer1_w0.646_ns1guess2.adj
les0822nh15t13y135_ocaer1_w0.646_ns1guess2.arad
les0822nh15t13y135_ocaer1_w0.646_ns1guess2.log
les0822nh15t13y135_ocaer1_w0.646_ns1true2.adj
les0822nh15t13y135_ocaer1_w0.646_ns1true2.arad
les0822nh15t13y135_ocaer1_w0.646_ns1true2.log
les0822nh15t13y135_w0.646_ns1guess2.adj
les0822nh15t13y135_w0.646_ns1guess2.arad
les0822nh15t13y135_w0.646_ns1guess2.log
les0822nh15t13y135_w0.646_ns1true2.adj
les0822nh15t13y135_w0.646_ns1true2.arad
les0822nh15t13y135_w0.646_ns1true2.log

Looking at the data files


In [4]:
# The log file
flog = "les0822nh15t13y135_w0.646_ns1true2.log"
print_head(flog, Nlines=10)

# The radiance file
frad = "les0822nh15t13y135_w0.646_ns1true2.arad"
print_head(frad, Nlines=25)

# The adjoint file
fadj_true = "les0822nh15t13y135_w0.646_ns1true2.adj"
fadj = "les0822nh15t13y135_w0.646_ns1guess2.adj"

    
print_head(fadj, Nlines=35)

# The adjoint Radiance file
farad = "les0822nh15t13y135_w0.646_ns1guess2.arad"
farad_true = "les0822nh15t13y135_w0.646_ns1true2.arad"
print_head(farad, Nlines=35)


Showing head: les0822nh15t13y135_w0.646_ns1true2.log
	
	Name of this SHDOMADJ run:
	les0822nh15t13y135_w0.646_ns1true2
	Input particle properties filename (or NONE for Rayleigh only)
	les0822nh15t13y135t.part
	Number of input scattering tables
	2
	Name of each scattering table
	water_w0.646.scat
	aero_oc_w0.646.scat


Showing head: les0822nh15t13y135_w0.646_ns1true2.arad
	! Polarized Spherical Harmonic Discrete Ordinate Radiative Transfer Radiance Output
	!  L= 15  M= 15  NLM=  256   NMU= 16  NPHI= 32  NANG=  354   NSH=   1767876
	!  NSTOKES=1   NX= 320   NY=   1   NZ=  83    NPTS=   36044   NCELLS=   42046
	!  PARTICLE_FILE=les0822nh15t13y135t.part
	!  CORRELATED_K-DIST_FILE=NONE   NUM_G= 1
	!  SOURCE_TYPE=SOLAR               DELTA-M METHOD
	!  GRID_TYPE=EVEN (X,Y)    PROPERTY-FILE (Z)    INDEPENDENT_PIXEL=2
	!  SURFACE_TYPE=FIXED LAMBERTIAN         HORIZ_BOUNDARY_COND=0
	!  GROUND_ALBEDO=0.0500000  SKY_RAD= 0.00000E+00
	!  SOLAR_FLUX= 0.100000E+01   SOLAR_MU=-0.9221000   SOLAR_AZ= 339.400
	!  UNITS=WATTS/(M^2 MICRON STER)    WAVELENGTH=      0.65
	!  SPLITTING_ACCURACY= 0.300E-01   SPHERICAL_HARMONIC_ACCURACY= 0.300E-02
	!  SOLUTION_ACCURACY= 0.300E-04
	!  MAXIMUM_ITERATIONS= 200   NUMBER_ITERATIONS= 123
	!  RADIANCE AT Z= 15.000
	!  320    1    1 ! NXOUT, NYOUT, NRADDIR
	!  ID      mu       phi
	!   1  1.0000000    0.00
	!
	ID  IY  IX    RADIANCE
	1   1   1   0.205754E+00
	1   1   2   0.196955E+00
	1   1   3   0.171690E+00
	1   1   4   0.196178E+00
	1   1   5   0.211501E+00


Showing head: les0822nh15t13y135_w0.646_ns1guess2.adj
	! SHDOMADJ optical property adjoint output
	!
	X       Y       Z      Extinct  SSalbedo  dMF/dExt  dMF/dSSalb  dMF/dP(Theta_j)
	0.0000  0.0000 -0.0087    0.00574 1.000000  0.000E+00  0.000E+00  0.000E+00
	0.0000  0.0000  0.0095    0.00574 1.000000  0.000E+00  0.000E+00  0.000E+00
	0.0000  0.0000  0.0315    0.00573 1.000000  0.000E+00  0.000E+00  0.000E+00
	0.0000  0.0000  0.0574    0.00571 1.000000  0.000E+00  0.000E+00  0.000E+00
	0.0000  0.0000  0.0878    0.00570 1.000000  0.000E+00  0.000E+00  0.000E+00
	0.0000  0.0000  0.1245    0.00568 1.000000  1.893E-12  0.000E+00  0.000E+00
	0.0000  0.0000  0.1650    0.00566 1.000000  9.660E-11  0.000E+00  0.000E+00
	0.0000  0.0000  0.2050    0.00564 1.000000  3.237E-10  0.000E+00  0.000E+00
	0.0000  0.0000  0.2450    0.00562 1.000000  5.613E-10  0.000E+00  0.000E+00
	0.0000  0.0000  0.2850    0.00561 1.000000  7.948E-10  0.000E+00  0.000E+00
	0.0000  0.0000  0.3250    0.00559 1.000000  8.362E-10  0.000E+00  0.000E+00
	0.0000  0.0000  0.3650    0.00557 1.000000  6.187E-10  0.000E+00  0.000E+00
	0.0000  0.0000  0.4050    0.00555 1.000000  3.825E-10  0.000E+00  0.000E+00
	0.0000  0.0000  0.4450    0.00553 1.000000  1.485E-10  0.000E+00  0.000E+00
	0.0000  0.0000  0.4850   11.03866 0.999998  7.969E-12 -6.166E-11  0.000E+00
	0.0000  0.0000  0.5250   21.63028 0.999998  2.652E-12 -4.458E-11  0.000E+00
	0.0000  0.0000  0.5650   30.82795 0.999997  3.485E-12 -8.711E-11  0.000E+00
	0.0000  0.0000  0.6050   38.37085 0.999997 -3.328E-10 -7.485E-08 -4.664E-06
	0.0000  0.0000  0.6450   40.02002 0.999997  3.989E-10 -3.205E-07 -1.682E-05
	0.0000  0.0000  0.6850   25.36786 0.999997  3.926E-09 -3.913E-07 -1.654E-05
	0.0000  0.0000  0.7250    0.00539 1.000000 -4.657E-08 -4.098E-10 -3.754E-09
	0.0000  0.0000  0.7650    0.00537 1.000000 -5.320E-09 -2.011E-10 -1.842E-09
	0.0000  0.0000  0.8050    0.00535 1.000000  1.510E-08 -7.357E-11 -6.740E-10
	0.0000  0.0000  0.8450    0.00532 1.000000  2.142E-08 -2.975E-11 -2.726E-10
	0.0000  0.0000  0.8850    0.00530 1.000000  2.389E-08 -1.395E-11 -1.278E-10
	0.0000  0.0000  0.9250    0.00528 1.000000  2.497E-08 -7.086E-12 -6.492E-11
	0.0000  0.0000  0.9650    6.33743 0.999998  1.209E-08 -4.835E-08 -3.338E-08
	0.0000  0.0000  1.0050   21.38868 0.999998  7.411E-09 -1.245E-07 -7.732E-08
	0.0000  0.0000  1.0450   36.03261 0.999997  1.077E-08 -3.283E-07 -4.714E-07
	0.0000  0.0000  1.0850   25.57642 0.999997  1.458E-08 -3.145E-07 -9.886E-07
	0.0000  0.0000  1.1250   14.30435 0.999998  1.149E-08 -1.339E-07 -6.860E-07
	0.0000  0.0000  1.1650   30.72805 0.999997  7.772E-09 -2.360E-07 -2.281E-06


Showing head: les0822nh15t13y135_w0.646_ns1guess2.arad
	! Polarized Spherical Harmonic Discrete Ordinate Radiative Transfer Radiance Output
	!  L= 15  M= 15  NLM=  256   NMU= 16  NPHI= 32  NANG=  354   NSH=   1684473
	!  NSTOKES=1   NX= 320   NY=   1   NZ=  83    NPTS=   35125   NCELLS=   40454
	!  PARTICLE_FILE=les0822nh15t13y135g.part
	!  CORRELATED_K-DIST_FILE=NONE   NUM_G= 1
	!  SOURCE_TYPE=SOLAR               DELTA-M METHOD
	!  GRID_TYPE=EVEN (X,Y)    PROPERTY-FILE (Z)    INDEPENDENT_PIXEL=2
	!  SURFACE_TYPE=FIXED LAMBERTIAN         HORIZ_BOUNDARY_COND=0
	!  GROUND_ALBEDO=0.0500000  SKY_RAD= 0.00000E+00
	!  SOLAR_FLUX= 0.100000E+01   SOLAR_MU=-0.9221000   SOLAR_AZ= 339.400
	!  UNITS=WATTS/(M^2 MICRON STER)    WAVELENGTH=      0.65
	!  SPLITTING_ACCURACY= 0.300E-01   SPHERICAL_HARMONIC_ACCURACY= 0.300E-02
	!  SOLUTION_ACCURACY= 0.300E-04
	!  MAXIMUM_ITERATIONS= 200   NUMBER_ITERATIONS= 123
	!  RADIANCE AT Z= 15.000
	!  320    1    1 ! NXOUT, NYOUT, NRADDIR
	!  ID      mu       phi
	!   1  1.0000000    0.00
	!
	ID  IY  IX    RADIANCE
	1   1   1   0.197190E+00
	1   1   2   0.188860E+00
	1   1   3   0.165215E+00
	1   1   4   0.187221E+00
	1   1   5   0.201468E+00
	1   1   6   0.204092E+00
	1   1   7   0.193736E+00
	1   1   8   0.170147E+00
	1   1   9   0.133220E+00
	1   1  10   0.107793E+00
	1   1  11   0.666059E-01
	1   1  12   0.487470E-01
	1   1  13   0.455724E-01
	1   1  14   0.446769E-01
	1   1  15   0.981936E-01



In [7]:
# Load the adjoint files into memory
adj_frame = pd.read_csv(fadj, quoting=3, delim_whitespace=True, skiprows=2)
adj_true_frame = pd.read_csv(fadj_true, quoting=3, delim_whitespace=True, skiprows=2)
adj_arad_frame = pd.read_csv(farad, quoting=3, delim_whitespace=True, skiprows=19)
adj_arad_true_frame = pd.read_csv(farad_true, quoting=3, delim_whitespace=True, skiprows=19)
    
    
# Get variables from the array
nx = adj_frame['X'].unique().size
nz = adj_frame["Z"].unique().size


# Get the guess case adjoint field properties
adj_x = np.array(adj_frame["X"], dtype='f8').reshape(nx, nz)
adj_z = np.array(adj_frame["Z"], dtype='f8').reshape(nx, nz)
adj_ext = np.array(adj_frame["Extinct"], dtype='f8').reshape(nx, nz)
adj_alb = np.array(adj_frame["SSalbedo"], dtype='f8').reshape(nx, nz)
adj_dMFdext = np.array(adj_frame["dMF/dExt"], dtype='f8').reshape(nx, nz)
adj_dMFdalb = np.array(adj_frame["dMF/dSSalb"], dtype='f8').reshape(nx, nz)

# Get the True case adjoint extinction field
adj_ext_true = np.array(adj_true_frame["Extinct"], dtype='f8').reshape(nx, nz) 

# Define the adjoint radiance parameters
adj_arad = np.array(adj_arad_frame['RADIANCE'], dtype='f8').reshape(nx)
adj_arad_true = np.array(adj_arad_true_frame['RADIANCE'], dtype='f8').reshape(nx)

In [8]:
# Check the sign of dMF/dExt for Z near the top of the atmosphere
adj_frame[abs(adj_frame["Z"]-15) < .2]


Out[8]:
X Y Z Extinct SSalbedo dMF/dExt dMF/dSSalb dMF/dP(Theta_j)
82 0.0000 0 15 0.00107 1 -5.7630 -0.007303 -0.066900
165 0.0625 0 15 0.00107 1 -5.4040 -0.006903 -0.063240
248 0.1250 0 15 0.00107 1 -4.5590 -0.005521 -0.050580
331 0.1875 0 15 0.00107 1 -6.1260 -0.007638 -0.069970
414 0.2500 0 15 0.00107 1 -6.7710 -0.008555 -0.078370
497 0.3125 0 15 0.00107 1 -6.5720 -0.008433 -0.077250
580 0.3750 0 15 0.00107 1 -5.7130 -0.007453 -0.068280
663 0.4375 0 15 0.00107 1 -4.9960 -0.006468 -0.059250
746 0.5000 0 15 0.00107 1 -3.4130 -0.004334 -0.039710
829 0.5625 0 15 0.00107 1 -2.9250 -0.003888 -0.035620
912 0.6250 0 15 0.00107 1 -2.7620 -0.003659 -0.033520
995 0.6875 0 15 0.00107 1 -2.0560 -0.002816 -0.025800
1078 0.7500 0 15 0.00107 1 -1.9100 -0.002647 -0.024250
1161 0.8125 0 15 0.00107 1 -1.7820 -0.002525 -0.023140
1244 0.8750 0 15 0.00107 1 -5.3360 -0.006965 -0.063810
1327 0.9375 0 15 0.00107 1 -6.7150 -0.008715 -0.079840
1410 1.0000 0 15 0.00107 1 -7.2750 -0.009091 -0.083280
1493 1.0625 0 15 0.00107 1 -7.3250 -0.009104 -0.083400
1576 1.1250 0 15 0.00107 1 -6.3700 -0.008180 -0.074930
1659 1.1875 0 15 0.00107 1 -5.0780 -0.006464 -0.059220
1742 1.2500 0 15 0.00107 1 -4.0340 -0.005292 -0.048480
1825 1.3125 0 15 0.00107 1 -4.3370 -0.005791 -0.053050
1908 1.3750 0 15 0.00107 1 -5.9460 -0.007700 -0.070540
1991 1.4375 0 15 0.00107 1 -6.5550 -0.008448 -0.077390
2074 1.5000 0 15 0.00107 1 -4.7300 -0.006188 -0.056690
2157 1.5625 0 15 0.00107 1 -4.4940 -0.005667 -0.051910
2240 1.6250 0 15 0.00107 1 -5.1200 -0.006369 -0.058340
2323 1.6875 0 15 0.00107 1 -5.4600 -0.006815 -0.062430
2406 1.7500 0 15 0.00107 1 -5.8380 -0.007253 -0.066440
2489 1.8125 0 15 0.00107 1 -4.9770 -0.006240 -0.057160
... ... ... ... ... ... ... ... ...
24152 18.1250 0 15 0.00107 1 -2.2390 -0.003262 -0.029880
24235 18.1875 0 15 0.00107 1 -1.7660 -0.002518 -0.023070
24318 18.2500 0 15 0.00107 1 -2.1470 -0.003023 -0.027690
24401 18.3125 0 15 0.00107 1 -2.9900 -0.003977 -0.036430
24484 18.3750 0 15 0.00107 1 -3.0620 -0.003977 -0.036430
24567 18.4375 0 15 0.00107 1 -5.6620 -0.007265 -0.066550
24650 18.5000 0 15 0.00107 1 -7.4070 -0.009339 -0.085550
24733 18.5625 0 15 0.00107 1 -5.8450 -0.007641 -0.070000
24816 18.6250 0 15 0.00107 1 -5.0250 -0.006580 -0.060280
24899 18.6875 0 15 0.00107 1 -4.6920 -0.006206 -0.056860
24982 18.7500 0 15 0.00107 1 -4.6610 -0.006121 -0.056070
25065 18.8125 0 15 0.00107 1 -3.8760 -0.004881 -0.044720
25148 18.8750 0 15 0.00107 1 -1.4780 -0.002024 -0.018540
25231 18.9375 0 15 0.00107 1 -1.9090 -0.002371 -0.021720
25314 19.0000 0 15 0.00107 1 -3.4870 -0.004342 -0.039770
25397 19.0625 0 15 0.00107 1 -4.0000 -0.005040 -0.046170
25480 19.1250 0 15 0.00107 1 -1.5030 -0.001863 -0.017060
25563 19.1875 0 15 0.00107 1 -0.2800 -0.000526 -0.004820
25646 19.2500 0 15 0.00107 1 0.3028 0.000045 0.000415
25729 19.3125 0 15 0.00107 1 -0.3874 -0.000854 -0.007822
25812 19.3750 0 15 0.00107 1 -1.8610 -0.002716 -0.024880
25895 19.4375 0 15 0.00107 1 -4.9990 -0.006609 -0.060540
25978 19.5000 0 15 0.00107 1 -6.3950 -0.008303 -0.076060
26061 19.5625 0 15 0.00107 1 -7.0560 -0.009055 -0.082950
26144 19.6250 0 15 0.00107 1 -6.6320 -0.008443 -0.077340
26227 19.6875 0 15 0.00107 1 -6.8060 -0.008879 -0.081340
26310 19.7500 0 15 0.00107 1 -6.8730 -0.008975 -0.082220
26393 19.8125 0 15 0.00107 1 -5.3910 -0.007234 -0.066270
26476 19.8750 0 15 0.00107 1 -5.8250 -0.007588 -0.069510
26559 19.9375 0 15 0.00107 1 -5.8220 -0.007439 -0.068150

320 rows × 8 columns


In [ ]:


In [9]:
# Plot uniformity parameters
H = 20 #5.5   # top of the atmosphere plots
mask_clouds = (adj_ext>= .75e-5) #* (adj_z <= 2.83)
cloud_outline = (adj_ext>= .75e-1) #* (adj_z <= 2.83)


# Scaling of the verticle coordinate 

Zscale = adj_z[0,:]
DZscale = 0 * Zscale
DZscale_other = 0 * Zscale
DZscale[:-1] += 1 / (Zscale[1:] - Zscale[:-1])
DZscale_other[1:] += 1 / (Zscale[1:] - Zscale[:-1])

alpha = .99725
DZscale = alpha * DZscale + (1-alpha) * DZscale_other



# Compute the residual
adj_residual = adj_arad_true - adj_arad

In [10]:
# Make a figure to plot the true extinction field
f0 = plt.figure(0, (15,8), facecolor='white')
ax0 = f0.add_axes([.1, .1, .8, .8])
contour0 = ax0.contourf(adj_x, adj_z, adj_ext_true, 100, cmap=cm.gray)
ax0.set_ybound((0,8))
ax0.set_title(" Guess for the extinction field at y=8.35")
f0.colorbar(contour0, ax=ax0)


Out[10]:
<matplotlib.colorbar.Colorbar instance at 0x4a63ab8>

In [11]:
# Make a figure to plot the true extinction field
f0 = plt.figure(0, (15,8), facecolor='white')
ax0 = f0.add_axes([.1, .1, .8, .8])

fvals = adj_ext_true - adj_ext

contour0 = ax0.contourf(adj_x, adj_z, fvals, 100, cmap=cm.gray)
ax0.set_ybound((0,8))
ax0.set_title("Difference $\sigma_{\mathrm{true}} - \sigma_{\mathrm{guess}}$")
f0.colorbar(contour0, ax=ax0)


Out[11]:
<matplotlib.colorbar.Colorbar instance at 0x5396dd0>

In [ ]:


In [12]:
# Make a figure to plot the true single-scattering albedo field
f1 = plt.figure(1, (15,8), facecolor='white')
ax1 = f1.add_axes([.1, .1, .8, .8])


levels = np.linspace(.999995, 1, 100)
contour1 = ax1.contourf(adj_x, adj_z, adj_alb, levels=levels, cmap=cm.gray)
ax1.set_ybound((0,13))
ax1.set_title("SSAlbedo field at y=8.35")
f1.colorbar(contour1, ax=ax1)


Out[12]:
<matplotlib.colorbar.Colorbar instance at 0x643f638>

In [13]:
# Make a figure to plot the true single-scattering albedo field
f2 = plt.figure(1, (10,6), facecolor='white')
ax2 = f2.add_axes([.1, .1, .8, .8])


fvals = adj_dMFdext * mask_clouds * DZscale * adj_ext# * .1 * adj_ext
fmax = abs(fvals).max() / 20.
levels = np.linspace(-fmax, fmax, 400)

contour2 = ax2.contourf(adj_x, adj_z, fvals, levels=levels, cmap=cm.RdBu_r, extend="both") #levels=linspace(-16, 16, 100),
contour2a = ax2.contour(adj_x, adj_z, cloud_outline, 1, cmap=cm.gray) #levels=linspace(-16, 16, 100),
#contour2b = ax2.contourf(adj_x, adj_z, cloud_outline * fvals, 10 , cmap=cm.RdBu_r, linewidth=3) #levels=linspace(-16, 16, 100),


# Add the residual to the plot
ftrans = lambda y: 200 * y + .8*H 
plot_overlay_origin = ftrans(0 * adj_residual)
plot_overlay_arad = ftrans(adj_residual)
plot_res_pos = ax2.fill_between(adj_x[:,0], plot_overlay_origin, plot_overlay_arad, 
                             where=(plot_overlay_origin<plot_overlay_arad), 
                             color='b')#, linewidth=3)
plot_res_neg = ax2.fill_between(adj_x[:,0], plot_overlay_origin, plot_overlay_arad, 
                             where=(plot_overlay_origin>plot_overlay_arad), 
                             color='r',)#, linewidth=3)
plot_res_line = ax2.plot(adj_x[:,0], plot_overlay_arad, 'k', linewidth=2)
plot_origin = ax2.plot(adj_x[:,0], plot_overlay_origin, 'k,', linewidth=1)
ax2.set_ybound((0,H))
ax2.set_ylabel("Z Height [km]", fontsize='x-large', fontweight='bold')
ax2.set_xlabel('X [km]', fontsize='large', fontweight='bold')
ax2.set_title("""
    Risidual ($y_{\mathrm{data}} - y_{\mathrm{model}}(\mathbf{a})$) and derivative (-$\mathrm{d}\Phi(\sigma) / \mathrm{d}\log\sigma$)
    """, fontsize='xx-large', fontweight='bold')
cb = f2.colorbar(contour2, ax=ax2)
cb.set_ticks([-fmax, 0, fmax])

f2.savefig("adjoint-derivative-single-scattering.jpg", dpi=300)

#cb.set_ticklabels(['-2', '-1', '0', '1', '2'])



In [14]:
# Make a figure to plot the true single-scattering albedo field
f2 = plt.figure(1, (10,6), facecolor='white')
ax2 = f2.add_axes([.1, .1, .8, .8])


fvals = -adj_dMFdext * mask_clouds * DZscale * adj_ext# * .1 * adj_ext

fvals = - (adj_ext_true - adj_ext) #* adj_ext

fmax = abs(fvals).max() / 1.
levels = np.linspace(-fmax, fmax, 200)

contour2 = ax2.contourf(adj_x, adj_z, fvals, levels=levels, cmap=cm.RdBu_r, extend="both") #levels=linspace(-16, 16, 100),
contour2a = ax2.contour(adj_x, adj_z, cloud_outline, 1, cmap=cm.gray) #levels=linspace(-16, 16, 100),
#contour2b = ax2.contourf(adj_x, adj_z, cloud_outline * fvals, 10 , cmap=cm.RdBu_r, linewidth=3) #levels=linspace(-16, 16, 100),


# Add the residual to the plot
ftrans = lambda y: 200 * y + .8*H 
plot_overlay_origin = ftrans(0 * adj_residual)
plot_overlay_arad = ftrans(adj_residual)
plot_res_pos = ax2.fill_between(adj_x[:,0], plot_overlay_origin, plot_overlay_arad, 
                             where=(plot_overlay_origin<plot_overlay_arad), 
                             color='b')#, linewidth=3)
plot_res_neg = ax2.fill_between(adj_x[:,0], plot_overlay_origin, plot_overlay_arad, 
                             where=(plot_overlay_origin>plot_overlay_arad), 
                             color='r',)#, linewidth=3)
plot_res_line = ax2.plot(adj_x[:,0], plot_overlay_arad, 'k', linewidth=2)
plot_origin = ax2.plot(adj_x[:,0], plot_overlay_origin, 'k,', linewidth=1)
ax2.set_ybound((0,H))
ax2.set_ylabel("Z Height [km]", fontsize='x-large', fontweight='bold')
ax2.set_xlabel('X [km]', fontsize='large', fontweight='bold')
ax2.set_title("""
    Risidual ($y_{\mathrm{data}} - y_{\mathrm{model}}(\mathbf{a})$) and extinction error ($\sigma_{\mathrm{true}}-\sigma_{\mathrm{guess}})$)
    """, fontsize='xx-large', fontweight='bold')
cb = f2.colorbar(contour2, ax=ax2)
cb.set_ticks([-fmax, 0, fmax])

f2.savefig("adjoint-extinction-difference-and-residual.jpg", dpi=300)

#cb.set_ticklabels(['-2', '-1', '0', '1', '2'])



In [15]:
plot(Zscale, DZscale)


Out[15]:
[<matplotlib.lines.Line2D at 0x406e050>]