<jupyter_start><jupyter_text>pyCRTM Jupyter Notebook  test_atms.ipynbDescription:Jupyter notebook version of the test_atms.ipynb Python test script for pyCRTM. Requirements:Python 3Record of Revisions:    Date: 2020-02-26 Author: <NAME> Description: Original code  Setup  Import Python utility modules:<jupyter_code>import configparser
import os
import h5py
import sys 
import numpy as np
from matplotlib import pyplot as plt<jupyter_output><empty_output><jupyter_text>Add current directory to Python sys path:<jupyter_code>thisDir = os.path.dirname(os.path.abspath('path'))
print(thisDir)
parentDir = os.path.dirname(thisDir)
print(parentDir)
sys.path.append(parentDir)<jupyter_output><empty_output><jupyter_text>Load the pyCRTM module:<jupyter_code>from pyCRTM import pyCRTM, profilesCreate<jupyter_output><empty_output><jupyter_text>Creating the Test Cases:<jupyter_code>thisDir = os.path.dirname(os.path.abspath('path'))
print(thisDir)
cases = os.listdir( os.path.join(thisDir,'data') )
print(cases)
cases.sort()<jupyter_output><empty_output><jupyter_text>Create 4 arrays to store atmospheric profile data and precomputed brightness temperatures:<jupyter_code>profiles = profilesCreate( 4, 92)
storedTb = []
storedEmis = []<jupyter_output><empty_output><jupyter_text>Populate the arrays with data for the 4 cases:<jupyter_code>for i,c in enumerate(cases):
        h5 = h5py.File(os.path.join(thisDir,'data',c) , 'r')
        profiles.Angles[i,0] = h5['zenithAngle'][()]
        profiles.Angles[i,1] = 999.9 
        profiles.Angles[i,2] = 100.0  # 100 degrees zenith below horizon.
        profiles.Angles[i,3] = 0.0 # zero solar azimuth 
        profiles.Angles[i,4] = h5['scanAngle'][()]
        profiles.DateTimes[i,0] = 2001
        profiles.DateTimes[i,1] = 1
        profiles.DateTimes[i,2] = 1
        profiles.Pi[i,:] = np.asarray(h5['pressureLevels'] )  # [hPa]
        profiles.P[i,:] = np.asarray(h5['pressureLayers'][()]) # [hPa]
        profiles.T[i,:] = np.asarray(h5['temperatureLayers']) # [K]
        profiles.Q[i,:] = np.asarray(h5['humidityLayers'])
        profiles.O3[i,:] = np.asarray(h5['ozoneConcLayers'])
        profiles.clouds[i,:,0,0] = np.asarray(h5['cloudConcentration'])
        profiles.clouds[i,:,0,1] = np.asarray(h5['cloudEffectiveRadius'])
        profiles.aerosols[i,:,0,0] = np.asarray(h5['aerosolConcentration'])
        profiles.aerosols[i,:,0,1] = np.asarray(h5['aerosolEffectiveRadius'])
        profiles.aerosolType[i] = h5['aerosolType'][()]
        profiles.cloudType[i] = h5['cloudType'][()]
        profiles.cloudFraction[i,:] = h5['cloudFraction'][()]
        profiles.climatology[i] = h5['climatology'][()]
        profiles.surfaceFractions[i,:] = h5['surfaceFractions']
        profiles.surfaceTemperatures[i,:] = h5['surfaceTemperatures']
        profiles.S2m[i,1] = 33.0 # just use salinity out of S2m for the moment.
        profiles.windSpeed10m[i] = 5.0
        profiles.LAI[i] = h5['LAI'][()]
        profiles.windDirection10m[i] = h5['windDirection10m'][()]
        # land, soil, veg, water, snow, ice
        profiles.surfaceTypes[i,0] = h5['landType'][()]
        profiles.surfaceTypes[i,1] = h5['soilType'][()]
        profiles.surfaceTypes[i,2] = h5['vegType'][()]
        profiles.surfaceTypes[i,3] = h5['waterType'][()]
        profiles.surfaceTypes[i,4] = h5['snowType'][()]
        profiles.surfaceTypes[i,5] = h5['iceType'][()]
        storedTb.append(np.asarray(h5['Tb_atms'][0:22]))
        storedEmis.append(np.asarray(h5['emissivity_atms'][0:22]))
        h5.close()<jupyter_output><empty_output><jupyter_text>Inspect input data  Plot the temperature:<jupyter_code>for i,c in enumerate(cases):
    plt.semilogy(profiles.T[i,:],profiles.P[i,:],label='case %d'%(i+1))
plt.gca().invert_yaxis()
plt.xlabel('Temperature [K]')
plt.ylabel('Pressure [hPa]')
plt.legend()<jupyter_output><empty_output><jupyter_text>Plot the water vapor concentration:<jupyter_code>for i,c in enumerate(cases):
    plt.loglog(profiles.Q[i,:],profiles.P[i,:],label='case %d'%(i+1))
plt.gca().invert_yaxis()
plt.xlabel(r'Absolute Humidity $[g/m^3]$')
plt.ylabel('Pressure [hPa]')
plt.legend()<jupyter_output><empty_output><jupyter_text>CRTM Initialization  Read the pyCRTM configuration file:<jupyter_code>pathInfo = configparser.ConfigParser()
pathInfo.read( os.path.join(parentDir,'crtm.cfg') )<jupyter_output><empty_output><jupyter_text>Read the path info of the CRTM binary coefficient files:<jupyter_code>coefficientPath = pathInfo['CRTM']['coeffs_dir']<jupyter_output><empty_output><jupyter_text>Set the satellite instrument ID as defined in the CRTM documentation:<jupyter_code>sensor_id = 'atms_npp'<jupyter_output><empty_output><jupyter_text>Instantiate a pyCRTM object:<jupyter_code>crtmOb = pyCRTM()<jupyter_output><empty_output><jupyter_text>Set the data members of the pyCRTM object:<jupyter_code>crtmOb.profiles = profiles
crtmOb.coefficientPath = pathInfo['CRTM']['coeffs_dir'] # Location of the CRTM binary coefficients
crtmOb.sensor_id = sensor_id
crtmOb.nThreads = 4<jupyter_output><empty_output><jupyter_text>Load CRTM SpcCoeff.bin Spectral coefficient files for the selected instrument:<jupyter_code>crtmOb.loadInst()<jupyter_output><empty_output><jupyter_text>Run the CRTM in Forward and Jacobian Mode Run the Forward solver of the CRTM for the 4 test cases:<jupyter_code>crtmOb.runDirect()<jupyter_output><empty_output><jupyter_text>Extract instrument brightness temperatures from solution object:<jupyter_code>forwardTb = crtmOb.Bt<jupyter_output><empty_output><jupyter_text>Extract surface emissivity from solution object:<jupyter_code>forwardEmissivity = crtmOb.surfEmisRefl[0,:]<jupyter_output><empty_output><jupyter_text>Reset the surface emissivity for the K-Matrix run:<jupyter_code>crtmOb.surfEmisRefl = []<jupyter_output><empty_output><jupyter_text>Run the K-Matrix solver:<jupyter_code>crtmOb.runK()<jupyter_output><empty_output><jupyter_text>Extract the brightness temperature solution and the surface emissivity Jacobian:<jupyter_code>kTb = crtmOb.Bt
kEmissivity = crtmOb.surfEmisRefl[0,:]<jupyter_output><empty_output><jupyter_text>Compare the extracted solution data with the precomputed Fortran solution from the CRTM and plot the difference if substantial:<jupyter_code>if ( all( np.abs( forwardTb.flatten() - np.asarray(storedTb).flatten() ) <= 1e-5)  and all( np.abs( kTb.flatten() - np.asarray(storedTb).flatten() ) <= 1e-5) ):
    print("Yay! all values are close enough to what CRTM test program produced!")
else: 
    print("Boo! something failed. Look at atms plots")
    wavenumbers = np.zeros([4,22])
    wavenumbers[0:4,:] = np.linspace(1,23,22)
        
    plt.figure()
    plt.plot(wavenumbers.T,forwardTb.T-np.asarray(storedTb).T ) 
    plt.legend(['1','2','3','4'])
    plt.savefig(os.path.join(thisDir,'atms'+'_spectrum_forward.png'))
    plt.figure()
    plt.plot(wavenumbers.T,forwardEmissivity.T-np.asarray(storedEmis).T)
    plt.savefig(os.path.join(thisDir,'atms'+'_emissivity_forward.png')) 
    
    plt.figure()
    plt.plot(wavenumbers.T,kTb.T-np.asarray(storedTb).T)
    plt.savefig(os.path.join(thisDir,'atms'+'_spectrum_k.png'))
    plt.figure()
    plt.plot(wavenumbers.T,kEmissivity.T-np.asarray(storedEmis).T)
    plt.savefig(os.path.join(thisDir,'atms'+'_emissivity_k.png')) 
    sys.exit("Boo! didn't pass tolerance with CRTM test program.")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14836875592256668, 'fraction_numerical': 0.017869229727900365, 'mean_word_length': 5.734731084776664, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 86, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '99643', 'license': 'Apache-2.0', 'n_tokens_mistral': 2863, 'n_tokens_neox': 2815, 'n_words': 427, 'path': 'testCases/test_atms.ipynb', 'repo_name': 'ejsimon9/pycrtm', 'size': '7392'}