"""
Python interface to NASA Earth-GRAM2010 software.
:NOTE:
You must manually set the 'exec' variable in this module.
:REQUIREMENTS:
NASA Earth-GRAM2010 software (available on request from NASA)
:doc:`atpy`
:doc:`numpy`
:TO-DO_LIST:
Automatically search path for executable, rather than require manual setting.
Be much smarter about scraping data from the ASCII files (e.g., BeautifulSoup?)
:HISTORY:
2011-10-19 13:10 IJMC: Created.
"""
# Identify the GRAM executable file:
exc = '/Users/ianc/proj/atmo/earthgram2010/GRAM_code/gram_E10.x'
# Load modules:
## TBD
import os
import pdb
import numpy as np
__version__ = '0.1'
[docs]class box:
"""Empty object container.
"""
# 2011-10-19 13:28 IJMC: Added to gram.py (from my spitzer.py)
def __init__(self):
return
def __class__(self):
return 'baseObject'
[docs]def findBreaks(lines, delim='-----'):
"""Find section breaks in a GRAM output text file.
:INPUTS:
lines : sequence of strings
Text lines; data contents to be parsed.
:OPTIONS:
delim : string
String which, if found within any string of 'lines',
indicates a new section.
:OUTPUT:
locs : sequence of ints
indices of each section break within lines.
:NOTES:
Note that no end-of-file index will be returned, unless the
proper delimiter is found there.
"""
# 2011-10-19 13:24 IJMC: Created
locs = []
for ii, line in enumerate(lines):
if line.find(delim) > -1:
locs.append(ii)
return locs
[docs]def readSpecies(filename, verbose=False):
"""Parse a GRAM2010 'species.txt' output file.
:INPUTS:
filename : str
filename of desired speciest.txt output file.
:OPTIONS:
verbose : bool
if True, display additional debugging information.
:OUTPUTS:
data structure
:NOTES:
Refer to NASA/TM-2011-216467, Appendix D.2, for the appropriate
formatting of the species.txt file.
Tested for GRAM v2.0 (May 2011); changes in data formatting may
cause some aspects of this to function improperly.
:SEE_ALSO:
:func:`readOutput`
"""
# 2011-10-19 13:15 IJMC: Created
def parseHeader(lines):
"""Parse species.txt output header lines, and return a data structure."""
# Strip blank lines:
newlines = []
for line in lines:
if len(line.strip())>0: newlines.append(line)
hdr = box()
try:
hdr.software = newlines[0].strip().strip('*').strip()
except:
hdr.software = 'Could not parse software from header text'
try:
hdr.version = newlines[1][newlines[1].lower().find('v'):newlines[1].find(',')]
except:
hdr.version = 'Could not parse version from header text'
try:
utc_date, utc_time, jd = \
[field[0:12].strip() for field in newlines[3].split('=')[1::]]
hdr.utc_date = utc_date
hdr.utc_time = utc_time
hdr.jd = float(jd)
except:
hdr.utc_date = 'Could not parse date from header text'
hdr.utc_time = 'Could not parse time from header text'
hdr.jd = 'Could not parse Julian Date from header text'
try:
f10_7, mf10_7, apInd = [field[0:12].strip() for field in newlines[4].split('=')[1::]]
hdr.f10_7 = float(f10_7)
hdr.mean_f10_7 = float(mf10_7)
hdr.ap_Index = float(apInd)
except:
hdr.f10_7 = 'Could not parse f10.7 from header text'
hdr.mean_f10_7 = 'Could not parse mean f10.7 from header text'
hdr.ap_Index = 'Could not parse ap Index from header text'
return hdr
def parseData(lines, dat=None, ind=None, nlev=None):
"""Parse a 7-line species.txt data section and add to dat."""
# To permit exec calls:
localBox = box
localnp = np
if dat is None:
dat = box()
species = []
species_pairs = []
for line in lines:
species_pairs.append([el.strip() for el in line[-11::].split('|')])
species = species + [el.strip() for el in line[-11::].split('|')]
for s in species:
exec 'dat.%s = localBox()' % s in locals()
exec 'dat.%s.concentration_ppmv = localnp.zeros(nlev, dtype=float)' % s in locals()
exec 'dat.%s.number_density_per_m3 = localnp.zeros(nlev, dtype=float)' % s in locals()
dat.species_pairs = species_pairs
if hasattr(dat, 'Tot'):
dat.Tot.MW = np.zeros(nlev, dtype=float)
# Define fields:
fields = ['height_km', 'lat_geocentric_deg', 'long_deg', \
'radius_km', 'lat_geodetic_deg', 'time_sec']
for f in fields:
exec 'dat.%s = localnp.zeros(nlev, dtype=float)' % f in locals()
# Get the data
# format is: s8s7s8ss9s9s|s9s9ss4s|s4
pos = [[1,9], [10,17], [18,26], [28,37], [38,47], [50,59], [60,69], [71,75], [76,80]]
for ii in range(7):
lvals = [lines[ii][p1:p2] for p1, p2 in pos]
if ii==0:
dat.height_km[ind] = ( float(lvals[0]))
dat.lat_geocentric_deg[ind] = (float(lvals[1]))
dat.long_deg[ind] = ( float(lvals[2]))
elif ii==1:
dat.radius_km[ind] = ( lvals[0])
dat.lat_geodetic_deg[ind] = (lvals[1])
elif ii==2:
dat.time_sec[ind] = (lvals[0])
if ii==6:
dat.N.concentration_ppmv[ind] = (float(lvals[3]))
dat.N.number_density_per_m3[ind] = (float(lvals[4]))
dat.Tot.MW[ind] = (float(lvals[5].strip('MW=')))
dat.Tot.number_density_per_m3[ind] = (float(lvals[6].strip('MW=')))
else:
exec 'dat.%s.concentration_ppmv[ind] = (float(lvals[3]))' % \
dat.species_pairs[ii][0] in locals()
exec 'dat.%s.number_density_per_m3[ind] = (float(lvals[4]))' % \
dat.species_pairs[ii][0] in locals()
exec 'dat.%s.concentration_ppmv[ind] = (float(lvals[5]))' % \
dat.species_pairs[ii][1] in locals()
exec 'dat.%s.number_density_per_m3[ind] = (float(lvals[6]))' % \
dat.species_pairs[ii][1] in locals()
return dat
out = box()
if os.path.isfile(filename):
f = open(filename, 'r')
raw = f.readlines()
f.close()
else:
print "File %s not found for opening by readSpecies."
breakLocs = findBreaks(raw, delim='----+----')
breakLocs.append(breakLocs[-1] + 8)
n_levels = len(breakLocs)
out.header = parseHeader(raw[0:breakLocs[0]])
out.data = None
# Determine how many levels there are:
n_lev = 0
for ii in range(len(breakLocs) - 1):
if (breakLocs[ii+1] - breakLocs[ii])==8:
n_lev += 1
out.header.n_lev = n_lev
for ii in range(n_lev):
if (breakLocs[ii+1] - breakLocs[ii])==8:
out.data = parseData(raw[breakLocs[ii]+1:breakLocs[ii+1]], \
dat=out.data, ind=ii, nlev=n_lev)
#
#out.dat = None
#for ii in range(len(breakLocs) - 1):
# out.dat = parseData(raw[breakLocs[ii]+1:breakLocs[ii+1]], dat=out.dat)
return out
[docs]def readOutput(filename, verbose=False):
"""Parse a GRAM2010 'output.txt' output file.
:INPUTS:
filename : str
filename of desired output.txt output file.
:OPTIONS:
verbose : bool
if True, display additional debugging information.
:OUTPUTS:
data structure
:NOTES:
Refer to NASA/TM-2011-216467, Appendix D.1, for the appropriate
formatting of the species.txt file.
Tested for GRAM v2.0 (May 2011); changes in data formatting may
cause some aspects of this to function improperly.
:TO_DO:
Much additional data: RRA sites, winds, variability parameters, etc.
:SEE_ALSO:
:func:`readSpecies`
"""
# 2011-10-19 20:21 IJMC: Created, from readSpecies
def parseHeader(lines):
"""Parse species.txt output header lines, and return a data structure."""
# Strip blank lines:
newlines = []
for line in lines:
if len(line.strip())>0: newlines.append(line)
hdr = box()
try:
hdr.software = newlines[0].strip().strip('*').strip()
except:
hdr.software = 'Could not parse software from header text'
try:
hdr.version = newlines[1][newlines[1].lower().find('v'):newlines[1].find(',')]
except:
hdr.version = 'Could not parse version from header text'
try:
utc_date, utc_time, jd = \
[field[0:12].strip() for field in newlines[2].split('=')[1::]]
hdr.utc_date = utc_date
hdr.utc_time = utc_time
hdr.jd = float(jd)
except:
hdr.utc_date = 'Could not parse date from header text'
hdr.utc_time = 'Could not parse time from header text'
hdr.jd = 'Could not parse Julian Date from header text'
try:
f10_7, mf10_7, apInd = [field[0:12].strip() for field in newlines[3].split('=')[1::]]
hdr.f10_7 = float(f10_7)
hdr.mean_f10_7 = float(mf10_7)
hdr.ap_Index = float(apInd)
except:
hdr.f10_7 = 'Could not parse f10.7 from header text'
hdr.mean_f10_7 = 'Could not parse mean f10.7 from header text'
hdr.ap_Index = 'Could not parse ap Index from header text'
return hdr
def parseData(lines, dat=None, ind=None, nlev=None):
"""Parse a 12-line species.txt data section and add to dat."""
# to permit calls to exec:
localBox = box
localnp = np
if dat is None:
dat = box()
dat.H2O = box()
dat.H2O.density_kg_m3 = np.zeros(nlev, dtype=float)
dat.H2O.pressure_pa = np.zeros(nlev, dtype=float)
# Define fields:
fields = ['height_km', 'lat_geocentric_deg', 'long_deg', \
'radius_km', 'lat_geodetic_deg', 'time_sec', \
'pressure_pa', 'density_kg_m3', 'temperature_K', \
'dewpoint_K', 'relative_humidity']
for f in fields:
exec 'dat.%s = localnp.zeros(nlev, dtype=float)' % f in locals()
# Get the data
# format is: s8s7s8ss9s9s|s9s9ss4s|s4
pos = [[1,9], [10,17], [18,26], [27,36], [37,47], \
[48,54], [55,61], [62,68], [69,75], [76,80]]
# Use mean values for pressure, density, temperature, RH.
# Skip much data, for now
for ii in [0, 1, 2, 10]:
lvals = [lines[ii][p1:p2] for p1, p2 in pos]
if ii==0:
dat.height_km[ind] = ( float(lvals[0]))
dat.lat_geocentric_deg[ind] = (float(lvals[1]))
dat.long_deg[ind] = ( float(lvals[2]))
dat.pressure_pa[ind] = ( float(lvals[3]))
dat.density_kg_m3[ind] = ( float(lvals[4]))
dat.temperature_K[ind] = ( float(lvals[5]))
elif ii==1:
dat.radius_km[ind] = ( lvals[0])
dat.lat_geodetic_deg[ind] = (lvals[1])
elif ii==2:
dat.time_sec[ind] = (lvals[0])
elif ii==10:
dat.H2O.pressure_pa[ind] = ( float(lvals[3]))
dat.H2O.density_kg_m3[ind] = (float(lvals[4]))
dat.dewpoint_K[ind] = ( float(lvals[5]))
dat.relative_humidity[ind] = (float(lvals[8].replace('%', '')))
return dat
out = box()
if os.path.isfile(filename):
f = open(filename, 'r')
raw = f.readlines()
f.close()
else:
print "File %s not found for opening by readOutput."
breakLocs = findBreaks(raw, delim='- -')[1::]
breakLocs.append(breakLocs[-1] + 13)
n_levels = len(breakLocs)
out.header = parseHeader(raw[0:breakLocs[0]])
out.data = None
# Determine how many levels there are:
n_lev = 0
for ii in range(len(breakLocs) - 1):
if (breakLocs[ii+1] - breakLocs[ii])==13:
n_lev += 1
out.header.n_lev = n_lev
for ii in range(n_lev):
if (breakLocs[ii+1] - breakLocs[ii])==13:
out.data = parseData(raw[breakLocs[ii]+1:breakLocs[ii+1]], \
dat=out.data, ind=ii, nlev=n_lev)
return out
[docs]def genAtm(outputFile, speciesFile, filename=None, \
clobber=False, species=['H2O', 'O3', 'N2O', 'CO', 'CH4']):
"""Generate .atm-format ASCII text from GRAM atmospheric data.
:INPUTS:
outputFile : str or box
GRAM 'output.txt' file name, or output of :func:`readOutput`.
If the latter, elements must have the same size as
speciesFile.
speciesFile : str or box
GRAM 'species.txt' file name, or output of :func:`readSpecies`.
If the latter, elements must have the same size as
outputFile.
:OPTIONS:
filename : str
Filename to write data to (if any)
clobber : bool
if True, overwrite existing filename (if any)
species : list of strings
species to write into output file.
:OUTPUTS:
lines : list of strings
Strings for writing to disk (or written to disk, if filename is not None)
:NOTES:
Formatting is as follows (note that some units are different than for GRAM):
%ALT PRES TEMP DENSITY H20 O3 N20 CO CH4
%(KM) (MB) (K) (CM-3) (PPMV) (PPMV) (PPMV) (PPMV) (PPMV)
0.00 1.013E+03 299.7 2.450E+19 2.59E+04 2.87E-02 3.20E-01 1.50E-01 1.70E+00
120.00 2.250E-05 380.0 4.225E+11 2.00E-01 5.00E-04 1.85E-04 5.00E+01 3.00E-02
I (IJMC) use these .atm files with D. Feldman's LBLRTM Matlab wrapper scripts
:SEE_ALSO:
:func:`readOutput`, :func:`readSpecies`
"""
# 2011-10-19 21:23 IJMC: Created
# Test inputs:
if not hasattr(outputFile, 'data'):
outputFile = readOutput(outputFile)
if not hasattr(speciesFile, 'data'):
speciesFile = readSpecies(speciesFile)
nspecies = 0
extant_species = []
for s in species:
if hasattr(speciesFile.data, s):
extant_species.append(s)
nspecies += 1
alt_index = np.argsort(outputFile.data.height_km)
# Create header lines:
lines = []
lines.append('%% Model atmosphere file from %s, %s\n' % \
(outputFile.header.software, outputFile.header.version) )
lines.append('%ALT PRES TEMP DENSITY' + \
''.join(['%7s ' % s for s in extant_species]) + '\n')
lines.append('%(KM) (MB) (K) (CM-3)' + (' (PPMV)'*nspecies) + '\n')
formatstr = '%6.2f %5.3e %5.1f %5.3e' + (' %4.2e' * nspecies) + '\n'
for ii in alt_index:
vals = (outputFile.data.height_km[ii], \
outputFile.data.pressure_pa[ii]/100., \
outputFile.data.temperature_K[ii], \
speciesFile.data.Tot.number_density_per_m3[ii]/1e6)
for s in extant_species:
vals += (eval('speciesFile.data.%s.concentration_ppmv[ii]' % s) ,)
lines.append(formatstr % vals)
# Write data:
if filename is not None:
if os.path.isfile(filename) and not clobber:
print "Could not write file '%s' to disk, because file already exists (try setting clobber)" % filename
else:
f = open(filename, 'w')
f.writelines(lines)
f.close()
return lines