sciapy.level2.density.main() - Code Metrics - st-bender/sciapy - Measure and Improve Code Quality continuously with Scrutinizer

sciapy.level2.density.main() A
last analyzed 2023-01-30 12:22 UTC

↳ Parent: sciapy.level2.density

Complexity

Conditions

Size

Total Lines	18
Code Lines	17

Duplication

Lines	0
Ratio	0 %

Code Coverage

Tests	1
CRAP Score	17.1835

Importance

Changes

Metric	Value
cc	4
eloc	17
nop	1
dl	0
loc	18
ccs	1
cts	16
cp	0.0625
crap	17.1835
rs	9.55
c	0
b	0
f	0

# -*- coding: utf-8 -*-
# vim:fileencoding=utf-8
#
# Copyright (c) 2015-2018 Stefan Bender
#
# This file is part of sciapy.
# sciapy is free software: you can redistribute it or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2.
# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html.
"""SCIAMACHY level 2 number density retrieval results interface

Interface classes for the level 2 retrieval results from text (ascii)
files and netcdf files for further processing.
"""

from __future__ import absolute_import, division, print_function

import os
import re
import sys
import datetime as dt

import numpy as np
try:
	from netCDF4 import Dataset as netcdf_file
	fmtargs = {"format": "NETCDF4"}
except ImportError:
	try:
		from scipy.io.netcdf import netcdf_file
		fmtargs = {"version": 1}
	except ImportError:
		from pupynere import netcdf_file
		fmtargs = {"version": 1}

__all__ = ["scia_densities", "_UTC"]

try:
	_UTC = dt.timezone.utc
except AttributeError:
	# python 2.7
	class UTC(dt.tzinfo):

		def utcoffset(self, d):
			return dt.timedelta(0)
		def tzname(self, d):
			return "UTC"
		def dst(self, d):
			return dt.timedelta(0)
	_UTC = UTC()

_meas_dtypes = [
	# initial output << v1.0
	[('gp_id', int),
		('alt_max', float), ('alt', float), ('alt_min', float),
		('lat_max', float), ('lat', float), ('lat_min', float),
		('density', float), ('dens_err_meas', float),
		('dens_err_tot', float), ('dens_tot', float)],
	# < 1.0 (NO_emiss-178-g729efb0)
	[('gp_id', int),
		('alt_max', float), ('alt', float), ('alt_min', float),
		('lat_max', float), ('lat', float), ('lat_min', float),
		('longitude', float),
		('density', float), ('dens_err_meas', float),
		('dens_err_tot', float), ('dens_tot', float)],
	# < 1.5 (NO_emiss-183-gcaa9349)
	[('gp_id', int),
		('alt_max', float), ('alt', float), ('alt_min', float),
		('lat_max', float), ('lat', float), ('lat_min', float),
		('longitude', float),
		('density', float), ('dens_err_meas', float),
		('dens_err_tot', float), ('dens_tot', float),
		('apriori', float)],
	# >= 1.5 (NO-v1.5)
	[('gp_id', int),
		('alt_max', float), ('alt', float), ('alt_min', float),
		('lat_max', float), ('lat', float), ('lat_min', float),
		('longitude', float),
		('density', float), ('dens_err_meas', float),
		('dens_err_tot', float), ('dens_tot', float),
		('apriori', float), ('akdiag', float)],
]


class scia_densities(object):
	"""SCIAMACHY orbital retrieved number densities

	Class interface to orbit-wise SCIAMACHY retrieval results.
	The attributes are based on the text file layout and are
	tied to the NO retrieval for now.

	Parameters
	----------
	ref_date: str, optional
		The reference date on which to base the date calculations on.
		Default: "2000-01-01"
	ver: str, optional
		Explicit density version, used for exporting the data.
		Not used if set to `None`.
		Default: `None`
	data_ver: str, optional
		Level 2 data version to use, as "ver" used for exporting.
		Not used if set to `None`.
		Default: `None`

	Attributes
	----------
	version
		file version string
	data_version
		level 2 data version
	date0
		reference date
	nalt
		number of altitudes in the orbit
	nlat
		number of latitudes in the orbit
	nlon
		number of longitudes in the orbit, if longitudes are available
	orbit
		SCIAMACHY/Envisat orbit number
	date
		number of days of the orbit counting from the reference date
		date0
	alts_min
	alts
	alts_max
		the altitude bins: minimum, central, and maximum altitude
	lats_min
	lats
	lats_max
		the latitude bins: minimum, central, and maximum latitude
	lons:
		the central longitude of the bins, only used if available

	densities
		NO number densities in the bins, (nlat, nalt) array_like
	dens_err_meas
		NO number densities measurement uncertainty,
		(nlat, nalt) array_like
	dens_err_tot
		NO number densities total uncertainty, (nlat, nalt) array_like
	dens_tot
		total number densities calculated and interpolated NRLMSIS-00
		values, (nlat, nalt) array_like

	apriori
		prior NO number densities, (nlat, nalt) array_like if available,
		otherwise `None`
	akdiag
		diagonal element of the averaging kernel matrix at the retrieval
		grid point. (nlat, nalt) array_like if available otherwise `None`

	Methods
	-------
	read_from_textfile
	read_from_netcdf
	read_from_file
	write_to_textfile
	write_to_netcdf

	Note
	----
	The variables are empty when initialized, use one of the
	`read_from_...()` methods to fill with actual data.
	"""
	def __init__(self, author="unknown", ref_date="2000-01-01", ver=None, data_ver=None):
		self.author = author
		self.version = ver
		self.data_version = data_ver
		self.date0 = dt.datetime.strptime(ref_date, "%Y-%m-%d").replace(tzinfo=_UTC)
		self.nalt = 0
		self.nlat = 0
		self.nlon = 0
		self.orbit = -1
		self.date = -1
		self.alts_min = np.array([])
		self.alts = np.array([])
		self.alts_max = np.array([])
		self.lats_min = np.array([])
		self.lats = np.array([])
		self.lats_max = np.array([])
		self.lons = np.array([])
		self.akdiag = None
		self.apriori = None

	def read_from_textfile(self, filename):
		"""Read NO densities from ascii table file

		Parameters
		----------
		filename: str, file object or io.TextIOBase.buffer
			The filename or stream to read the data from. For example
			to read from stdin in python 3, pass `sys.stdin.buffer`.
		"""
		def _unsrt_unique(a):
			return a[np.sort(np.unique(a, return_index=True)[1])]

		if hasattr(filename, 'seek'):
			f = filename
		else:
			f = open(filename, 'rb')
			# example filename:000NO_orbit_41467_20100203_Dichten.txt
			fn_fields = os.path.basename(filename).split('_')
			self.orbit = int(fn_fields[2])
			self.date = (dt.datetime.strptime(fn_fields[3], "%Y%m%d")
						.replace(tzinfo=_UTC) - self.date0).days
			if self.data_version is None:
				# try some heuristics to find the level 2 data version
				_dir = os.path.dirname(filename)
				_m = re.search(".*[_-]v([0-9]+[.].*)", _dir)
				if _m:
					self.data_version = _m.group(1)
				else:
					self.data_version = "unknown"
		data = f.readline().split()
		mydtype = _meas_dtypes[len(data) - 13]
		marr = np.genfromtxt(f, dtype=mydtype)
		f.close()

		# unique altitudes
		self.alts_min = _unsrt_unique(marr['alt_min'])
		self.alts = _unsrt_unique(marr['alt'])
		self.alts_max = _unsrt_unique(marr['alt_max'])

		# unique latitudes
		self.lats_min = _unsrt_unique(marr['lat_min'])
		self.lats = _unsrt_unique(marr['lat'])
		self.lats_max = _unsrt_unique(marr['lat_max'])

		# unique longitudes if available
		try:
			self.lons = _unsrt_unique(marr['longitude'])
			self.nlon = len(self.lons)
		except ValueError:
			pass

		self.nalt = len(self.alts)
		self.nlat = len(self.lats)

		# reorder by latitude first, then altitude
		self.densities = marr['density'].flatten().reshape(self.nalt, self.nlat).transpose()
		self.dens_err_meas = marr['dens_err_meas'].flatten().reshape(self.nalt, self.nlat).transpose()
		self.dens_err_tot = marr['dens_err_tot'].flatten().reshape(self.nalt, self.nlat).transpose()
		self.dens_tot = marr['dens_tot'].flatten().reshape(self.nalt, self.nlat).transpose()

		# apriori if available
		try:
			self.apriori = marr['apriori'].flatten().reshape(self.nalt, self.nlat).transpose()
		except ValueError:
			pass
		# akdiag if available
		try:
			self.akdiag = marr['akdiag'].flatten().reshape(self.nalt, self.nlat).transpose()
		except ValueError:
			pass

	def write_to_textfile(self, filename):
		"""Write NO densities to ascii table files

		Parameters
		----------
		filename: str or file object or io.TextIOBase.buffer
			The filename or stream to write the data to. For writing to
			stdout in python 3, pass `sys.stdout.buffer`.
		"""
		if hasattr(filename, 'seek'):
			f = filename
		else:
			f = open(filename, 'w')

		header = "%5s %13s %12s %13s %13s %12s %13s  %13s %12s %12s %12s" % ("GP_ID",
				"Max_Hoehe[km]", "Hoehe[km]", "Min_Hoehe[km]",
				"Max_Breite[°]", "Breite[°]", "Min_Breite[°]",
				"Dichte[cm^-3]", "Fehler Mess[cm^-3]",
				"Fehler tot[cm^-3]", "Gesamtdichte[cm^-3]")
		if self.nlon > 0:
			header = header[:87] + " %13s" % ("Laenge[°]",) + header[87:]
		if self.apriori is not None:
			header = header + " %12s" % ("apriori[cm^-3]",)
		if self.akdiag is not None:
			header = header + " %12s" % ("AKdiag",)
		print(header, file=f)

		oformat = "%5i  %+1.5E %+1.5E  %+1.5E  %+1.5E %+1.5E  %+1.5E   %+1.5E       %+1.5E      %+1.5E        %+1.5E"
		if self.nlon > 0:
			oformat = oformat[:49] + "  %+1.5E" + oformat[49:]
		oformata = "  %+1.5E"

		for i, a in enumerate(self.alts):
			for j, b in enumerate(self.lats):
				line_list = [i * self.nlat + j,
					self.alts_max[i], a, self.alts_min[i],
					self.lats_max[j], b, self.lats_min[j],
					self.densities[j, i], self.dens_err_meas[j, i],
					self.dens_err_tot[j, i], self.dens_tot[j, i]]
				if self.nlon > 0:
					line_list.insert(7, self.lons[j])
				print(oformat % tuple(line_list),
					end="", file=f)
				if self.apriori is not None:
					print(" " + oformata % self.apriori[j, i], end="", file=f)
				if self.akdiag is not None:
					print(" " + oformata % self.akdiag[j, i], end="", file=f)
				print("", file=f)

	def write_to_netcdf(self, filename, close=True):
		"""Write NO densities to netcdf files

		This function has no stream, i.e. file object, support.

		Parameters
		----------
		filename: str
			The name of the file to write the data to.
		close: bool, optional
			Whether or not to close the file after writing.
			Setting to `False` enables appending further data
			to the same file.
			Default: True

		Returns
		-------
		Nothing if `close` is True. If `close` is False, returns either an
		`netCDF4.Dataset`,
		`scipy.io.netcdf.netcdf_file` or
		`pupynere.netcdf_file` instance depending on availability.
		"""
		alts_min_out = np.asarray(self.alts_min).reshape(self.nalt)
		alts_out = np.asarray(self.alts).reshape(self.nalt)
		alts_max_out = np.asarray(self.alts_max).reshape(self.nalt)

		lats_min_out = np.asarray(self.lats_min).reshape(self.nlat)
		lats_out = np.asarray(self.lats).reshape(self.nlat)
		lats_max_out = np.asarray(self.lats_max).reshape(self.nlat)

		ncf = netcdf_file(filename, 'w', **fmtargs)

		if self.version is not None:
			ncf.version = self.version
		if self.data_version is not None:
			ncf.L2_data_version = self.data_version
		#ncf.creation_time = dt.datetime.utcnow().replace(tzinfo=_UTC).strftime("%a %b %d %Y %H:%M:%S %z (%Z)")
		ncf.creation_time = dt.datetime.utcnow().strftime("%a %b %d %Y %H:%M:%S +00:00 (UTC)")
		ncf.author = self.author

		# create netcdf file
		ncf.createDimension('time', None)
		ncf.createDimension('altitude', self.nalt)
		ncf.createDimension('latitude', self.nlat)

		forbit = ncf.createVariable('orbit', np.dtype('int32').char, ('time',))
		ftime = ncf.createVariable('time', np.dtype('int32').char, ('time',))

		falts_min = ncf.createVariable('alt_min', np.dtype('float64').char, ('altitude',))
		falts = ncf.createVariable('altitude', np.dtype('float64').char, ('altitude',))
		falts_max = ncf.createVariable('alt_max', np.dtype('float64').char, ('altitude',))
		flats_min = ncf.createVariable('lat_min', np.dtype('float64').char, ('latitude',))
		flats = ncf.createVariable('latitude', np.dtype('float64').char, ('latitude',))
		flats_max = ncf.createVariable('lat_max', np.dtype('float64').char, ('latitude',))

		falts_min.units = 'km'
		falts_min.positive = 'up'
		falts.units = 'km'
		falts.positive = 'up'
		falts_max.units = 'km'
		falts_max.positive = 'up'
		flats_min.units = 'degrees_north'
		flats.units = 'degrees_north'
		flats_max.units = 'degrees_north'

		forbit.units = '1'
		forbit.long_name = 'SCIAMACHY/Envisat orbit number'
		ftime.units = 'days since {0}'.format(self.date0.isoformat(sep=' '))
		ftime.standard_name = 'time'

		fdens = ncf.createVariable('density', np.dtype('float64').char, ('time', 'latitude', 'altitude'))
		fdens.units = 'cm^{-3}'
		fdens.standard_name = 'number_concentration_of_nitrogen_monoxide_molecules_in_air'
		fdens_err_meas = ncf.createVariable('error_meas', np.dtype('float64').char, ('time', 'latitude', 'altitude'))
		fdens_err_meas.units = 'cm^{-3}'
		fdens_err_meas.long_name = 'NO number density measurement error'
		fdens_err_tot = ncf.createVariable('error_tot', np.dtype('float64').char, ('time', 'latitude', 'altitude'))
		fdens_err_tot.units = 'cm^{-3}'
		fdens_err_tot.long_name = 'NO number density total error'
		fdens_tot = ncf.createVariable('density_air', np.dtype('float64').char, ('time', 'latitude', 'altitude'))
		fdens_tot.units = 'cm^{-3}'
		fdens_tot.long_name = 'approximate overall number concentration of air molecules (NRLMSIS-00)'

		ftime[:] = [self.date]
		forbit[:] = [self.orbit]

		falts_min[:] = alts_min_out
		falts[:] = alts_out
		falts_max[:] = alts_max_out
		flats_min[:] = lats_min_out
		flats[:] = lats_out
		flats_max[:] = lats_max_out
		# reorder by latitude first, then altitude
		fdens[0, :] = self.densities
		# reorder by latitude first, then altitude
		fdens_err_meas[0, :] = self.dens_err_meas
		fdens_err_tot[0, :] = self.dens_err_tot
		fdens_tot[0, :] = self.dens_tot

		# longitudes if they are available
		if self.nlon > 0:
			lons_out = np.asarray(self.lons).reshape(self.nlon)
			flons = ncf.createVariable('longitude', np.dtype('float64').char, ('time', 'latitude',))
			flons.units = 'degrees_east'
			flons[0, :] = lons_out

		if self.apriori is not None:
			fapriori = ncf.createVariable('apriori',
					np.dtype('float64').char, ('time', 'latitude', 'altitude'))
			fapriori.units = 'cm^{-3}'
			fapriori.long_name = 'apriori NO number density'
			fapriori[0, :] = self.apriori

		if self.akdiag is not None:
			fakdiag = ncf.createVariable('akm_diagonal',
					np.dtype('float64').char, ('time', 'latitude', 'altitude'))
			fakdiag.units = '1'
			fakdiag.long_name = 'averaging kernel matrix diagonal element'
			fakdiag[0, :] = self.akdiag
		if close:
			ncf.close()
		else:
			return ncf

	def read_from_netcdf(self, filename, close=True):
		"""Read NO densities from netcdf files

		This function has no stream, i.e. file object support.

		Parameters
		----------
		filename: str
			The filename to read the data from.
		close: bool, optional
			Whether or not to close the file after reading.
			Setting to `False` enables reading further data
			from the same file.
			Default: True

		Returns
		-------
		Nothing if `close` is True. If `close` is False, returns either an
		`netCDF4.Dataset`,
		`scipy.io.netcdf.netcdf_file` or
		`pupynere.netcdf_file` instance depending on availability.
		"""
		def _try_decode(s):
			if hasattr(s, "decode"):
				return s.decode()
			return s
		ncf = netcdf_file(filename, 'r')

		try:
			self.author = _try_decode(ncf.author)
		except AttributeError:
			pass
		try:
			self.version = _try_decode(ncf.version)
		except AttributeError:
			pass
		try:
			self.data_version = _try_decode(ncf.L2_data_version)
		except AttributeError:
			pass

		self.alts_min = ncf.variables['alt_min'][:].copy()
		self.alts = ncf.variables['altitude'][:].copy()
		self.alts_max = ncf.variables['alt_max'][:].copy()
		self.lats_min = ncf.variables['lat_min'][:].copy()
		self.lats = ncf.variables['latitude'][:].copy()
		self.lats_max = ncf.variables['lat_max'][:].copy()

		self.nalt = len(self.alts)
		self.nlat = len(self.lats)

		self.date = ncf.variables['time'][:].copy()
		self.orbit = ncf.variables['orbit'][:].copy()

		self.densities = ncf.variables['density'][:].copy()
		self.dens_err_meas = ncf.variables['error_meas'][:].copy()
		self.dens_err_tot = ncf.variables['error_tot'][:].copy()
		self.dens_tot = ncf.variables['density_air'][:].copy()

		# longitudes if they are available
		try:
			self.lons = ncf.variables['longitude'][:].copy()
			self.nlon = self.lons.shape[1]
		except KeyError:
			pass

		# apriori
		try:
			self.apriori = ncf.variables['apriori'][:].copy()
		except KeyError:
			pass

		# akm diagonal elements
		try:
			self.akdiag = ncf.variables['akm_diagonal'][:].copy()
		except KeyError:
			pass

		if close:
			ncf.close()
		else:
			return ncf

	def read_from_file(self, filename):
		"""Wrapper to read NO desnities from files

		Simple wrapper to delegate reading the data from either netcdf
		or ascii files. Poor man's logic: simply try netcdf first, and
		if that fails, read as ascii.

		Parameters
		----------
		filename: str
			The filename to read the data from.
		"""
		try:
			# try netcdf first
			self.read_from_netcdf(filename)
		except (IOError, OSError, TypeError):
			# fall back to text file as a last resort
			self.read_from_textfile(filename)


def main(*args):
	argc = len(sys.argv)
	if argc < 2:
		print("Not enough arguments, Usage:\n"
			"{0} [input] output [< input]".format(sys.argv[0]))
		sys.exit(1)
	elif argc < 3:
		try:
			infile = sys.stdin.buffer  # Python 3
		except AttributeError:
			infile = sys.stdin
		outfile = sys.argv[1]
	else:
		infile = sys.argv[1]
		outfile = sys.argv[2]
	sdl = scia_densities()
	sdl.read_from_file(infile)

	sdl.write_to_netcdf(outfile)



if __name__ == "__main__":
	sys.exit(main())


1		# -- coding: utf-8 --
2		# vim:fileencoding=utf-8
3		#
4		# Copyright (c) 2015-2018 Stefan Bender
5		#
6		# This file is part of sciapy.
7		# sciapy is free software: you can redistribute it or modify it
8		# under the terms of the GNU General Public License as published by
9		# the Free Software Foundation, version 2.
10		# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html.
11	1	"""SCIAMACHY level 2 number density retrieval results interface
12
13		Interface classes for the level 2 retrieval results from text (ascii)
14		files and netcdf files for further processing.
15		"""
16
17	1	from __future__ import absolute_import, division, print_function
18
19	1	import os
20	1	import re
21	1	import sys
22	1	import datetime as dt
23
24	1	import numpy as np
25	1	try:
26	1	from netCDF4 import Dataset as netcdf_file
27	1	fmtargs = {"format": "NETCDF4"}
28		except ImportError:
29		try:
30		from scipy.io.netcdf import netcdf_file
31		fmtargs = {"version": 1}
32		except ImportError:
33		from pupynere import netcdf_file
34		fmtargs = {"version": 1}
35
36	1	__all__ = ["scia_densities", "_UTC"]
37
38	1	try:
39	1	_UTC = dt.timezone.utc
40		except AttributeError:
41		# python 2.7
42		class UTC(dt.tzinfo):
		0 ignored issues – show Comprehensibility Best Practice introduced 2018-08-13 16:18 UTC by Report Bug Copy Issue Report The variable `dt` does not seem to be defined. Loading history...
43		def utcoffset(self, d):
44		return dt.timedelta(0)
45		def tzname(self, d):
46		return "UTC"
47		def dst(self, d):
48		return dt.timedelta(0)
49		_UTC = UTC()
50
51	1	_meas_dtypes = [
52		# initial output << v1.0
53		[('gp_id', int),
54		('alt_max', float), ('alt', float), ('alt_min', float),
55		('lat_max', float), ('lat', float), ('lat_min', float),
56		('density', float), ('dens_err_meas', float),
57		('dens_err_tot', float), ('dens_tot', float)],
58		# < 1.0 (NO_emiss-178-g729efb0)
59		[('gp_id', int),
60		('alt_max', float), ('alt', float), ('alt_min', float),
61		('lat_max', float), ('lat', float), ('lat_min', float),
62		('longitude', float),
63		('density', float), ('dens_err_meas', float),
64		('dens_err_tot', float), ('dens_tot', float)],
65		# < 1.5 (NO_emiss-183-gcaa9349)
66		[('gp_id', int),
67		('alt_max', float), ('alt', float), ('alt_min', float),
68		('lat_max', float), ('lat', float), ('lat_min', float),
69		('longitude', float),
70		('density', float), ('dens_err_meas', float),
71		('dens_err_tot', float), ('dens_tot', float),
72		('apriori', float)],
73		# >= 1.5 (NO-v1.5)
74		[('gp_id', int),
75		('alt_max', float), ('alt', float), ('alt_min', float),
76		('lat_max', float), ('lat', float), ('lat_min', float),
77		('longitude', float),
78		('density', float), ('dens_err_meas', float),
79		('dens_err_tot', float), ('dens_tot', float),
80		('apriori', float), ('akdiag', float)],
81		]
82
83
84	1	class scia_densities(object):
85		"""SCIAMACHY orbital retrieved number densities
86
87		Class interface to orbit-wise SCIAMACHY retrieval results.
88		The attributes are based on the text file layout and are
89		tied to the NO retrieval for now.
90
91		Parameters
92		----------
93		ref_date: str, optional
94		The reference date on which to base the date calculations on.
95		Default: "2000-01-01"
96		ver: str, optional
97		Explicit density version, used for exporting the data.
98		Not used if set to `None`.
99		Default: `None`
100		data_ver: str, optional
101		Level 2 data version to use, as "ver" used for exporting.
102		Not used if set to `None`.
103		Default: `None`
104
105		Attributes
106		----------
107		version
108		file version string
109		data_version
110		level 2 data version
111		date0
112		reference date
113		nalt
114		number of altitudes in the orbit
115		nlat
116		number of latitudes in the orbit
117		nlon
118		number of longitudes in the orbit, if longitudes are available
119		orbit
120		SCIAMACHY/Envisat orbit number
121		date
122		number of days of the orbit counting from the reference date
123		date0
124		alts_min
125		alts
126		alts_max
127		the altitude bins: minimum, central, and maximum altitude
128		lats_min
129		lats
130		lats_max
131		the latitude bins: minimum, central, and maximum latitude
132		lons:
133		the central longitude of the bins, only used if available
134
135		densities
136		NO number densities in the bins, (nlat, nalt) array_like
137		dens_err_meas
138		NO number densities measurement uncertainty,
139		(nlat, nalt) array_like
140		dens_err_tot
141		NO number densities total uncertainty, (nlat, nalt) array_like
142		dens_tot
143		total number densities calculated and interpolated NRLMSIS-00
144		values, (nlat, nalt) array_like
145
146		apriori
147		prior NO number densities, (nlat, nalt) array_like if available,
148		otherwise `None`
149		akdiag
150		diagonal element of the averaging kernel matrix at the retrieval
151		grid point. (nlat, nalt) array_like if available otherwise `None`
152
153		Methods
154		-------
155		read_from_textfile
156		read_from_netcdf
157		read_from_file
158		write_to_textfile
159		write_to_netcdf
160
161		Note
162		----
163		The variables are empty when initialized, use one of the
164		`read_from_...()` methods to fill with actual data.
165		"""
166	1	def __init__(self, author="unknown", ref_date="2000-01-01", ver=None, data_ver=None):
167	1	self.author = author
168	1	self.version = ver
169	1	self.data_version = data_ver
170	1	self.date0 = dt.datetime.strptime(ref_date, "%Y-%m-%d").replace(tzinfo=_UTC)
171	1	self.nalt = 0
172	1	self.nlat = 0
173	1	self.nlon = 0
174	1	self.orbit = -1
175	1	self.date = -1
176	1	self.alts_min = np.array([])
177	1	self.alts = np.array([])
178	1	self.alts_max = np.array([])
179	1	self.lats_min = np.array([])
180	1	self.lats = np.array([])
181	1	self.lats_max = np.array([])
182	1	self.lons = np.array([])
183	1	self.akdiag = None
184	1	self.apriori = None
185
186	1	def read_from_textfile(self, filename):
187		"""Read NO densities from ascii table file
188
189		Parameters
190		----------
191		filename: str, file object or io.TextIOBase.buffer
192		The filename or stream to read the data from. For example
193		to read from stdin in python 3, pass `sys.stdin.buffer`.
194		"""
195	1	def _unsrt_unique(a):
196	1	return a[np.sort(np.unique(a, return_index=True)[1])]
197
198	1	if hasattr(filename, 'seek'):
199		f = filename
200		else:
201	1	f = open(filename, 'rb')
202		# example filename:000NO_orbit_41467_20100203_Dichten.txt
203	1	fn_fields = os.path.basename(filename).split('_')
204	1	self.orbit = int(fn_fields[2])
205	1	self.date = (dt.datetime.strptime(fn_fields[3], "%Y%m%d")
206		.replace(tzinfo=_UTC) - self.date0).days
207	1	if self.data_version is None:
208		# try some heuristics to find the level 2 data version
209	1	_dir = os.path.dirname(filename)
210	1	_m = re.search(".[_-]v([0-9]+[.].)", _dir)
211	1	if _m:
212	1	self.data_version = _m.group(1)
213		else:
214	1	self.data_version = "unknown"
215	1	data = f.readline().split()
216	1	mydtype = _meas_dtypes[len(data) - 13]
217	1	marr = np.genfromtxt(f, dtype=mydtype)
218	1	f.close()
219
220		# unique altitudes
221	1	self.alts_min = _unsrt_unique(marr['alt_min'])
222	1	self.alts = _unsrt_unique(marr['alt'])
223	1	self.alts_max = _unsrt_unique(marr['alt_max'])
224
225		# unique latitudes
226	1	self.lats_min = _unsrt_unique(marr['lat_min'])
227	1	self.lats = _unsrt_unique(marr['lat'])
228	1	self.lats_max = _unsrt_unique(marr['lat_max'])
229
230		# unique longitudes if available
231	1	try:
232	1	self.lons = _unsrt_unique(marr['longitude'])
233	1	self.nlon = len(self.lons)
234	1	except ValueError:
235	1	pass
236
237	1	self.nalt = len(self.alts)
238	1	self.nlat = len(self.lats)
239
240		# reorder by latitude first, then altitude
241	1	self.densities = marr['density'].flatten().reshape(self.nalt, self.nlat).transpose()
242	1	self.dens_err_meas = marr['dens_err_meas'].flatten().reshape(self.nalt, self.nlat).transpose()
243	1	self.dens_err_tot = marr['dens_err_tot'].flatten().reshape(self.nalt, self.nlat).transpose()
244	1	self.dens_tot = marr['dens_tot'].flatten().reshape(self.nalt, self.nlat).transpose()
245
246		# apriori if available
247	1	try:
248	1	self.apriori = marr['apriori'].flatten().reshape(self.nalt, self.nlat).transpose()
249	1	except ValueError:
250	1	pass
251		# akdiag if available
252	1	try:
253	1	self.akdiag = marr['akdiag'].flatten().reshape(self.nalt, self.nlat).transpose()
254	1	except ValueError:
255	1	pass
256
257	1	def write_to_textfile(self, filename):
258		"""Write NO densities to ascii table files
259
260		Parameters
261		----------
262		filename: str or file object or io.TextIOBase.buffer
263		The filename or stream to write the data to. For writing to
264		stdout in python 3, pass `sys.stdout.buffer`.
265		"""
266	1	if hasattr(filename, 'seek'):
267		f = filename
268		else:
269	1	f = open(filename, 'w')
270
271	1	header = "%5s %13s %12s %13s %13s %12s %13s %13s %12s %12s %12s" % ("GP_ID",
272		"Max_Hoehe[km]", "Hoehe[km]", "Min_Hoehe[km]",
273		"Max_Breite[°]", "Breite[°]", "Min_Breite[°]",
274		"Dichte[cm^-3]", "Fehler Mess[cm^-3]",
275		"Fehler tot[cm^-3]", "Gesamtdichte[cm^-3]")
276	1	if self.nlon > 0:
277	1	header = header[:87] + " %13s" % ("Laenge[°]",) + header[87:]
278	1	if self.apriori is not None:
279	1	header = header + " %12s" % ("apriori[cm^-3]",)
280	1	if self.akdiag is not None:
281	1	header = header + " %12s" % ("AKdiag",)
282	1	print(header, file=f)
283
284	1	oformat = "%5i %+1.5E %+1.5E %+1.5E %+1.5E %+1.5E %+1.5E %+1.5E %+1.5E %+1.5E %+1.5E"
285	1	if self.nlon > 0:
286	1	oformat = oformat[:49] + " %+1.5E" + oformat[49:]
287	1	oformata = " %+1.5E"
288
289	1	for i, a in enumerate(self.alts):
290	1	for j, b in enumerate(self.lats):
291	1	line_list = [i * self.nlat + j,
292		self.alts_max[i], a, self.alts_min[i],
293		self.lats_max[j], b, self.lats_min[j],
294		self.densities[j, i], self.dens_err_meas[j, i],
295		self.dens_err_tot[j, i], self.dens_tot[j, i]]
296	1	if self.nlon > 0:
297	1	line_list.insert(7, self.lons[j])
298	1	print(oformat % tuple(line_list),
299		end="", file=f)
300	1	if self.apriori is not None:
301	1	print(" " + oformata % self.apriori[j, i], end="", file=f)
302	1	if self.akdiag is not None:
303	1	print(" " + oformata % self.akdiag[j, i], end="", file=f)
304	1	print("", file=f)
305
306	1	def write_to_netcdf(self, filename, close=True):
307		"""Write NO densities to netcdf files
308
309		This function has no stream, i.e. file object, support.
310
311		Parameters
312		----------
313		filename: str
314		The name of the file to write the data to.
315		close: bool, optional
316		Whether or not to close the file after writing.
317		Setting to `False` enables appending further data
318		to the same file.
319		Default: True
320
321		Returns
322		-------
323		Nothing if `close` is True. If `close` is False, returns either an
324		`netCDF4.Dataset`,
325		`scipy.io.netcdf.netcdf_file` or
326		`pupynere.netcdf_file` instance depending on availability.
327		"""
328	1	alts_min_out = np.asarray(self.alts_min).reshape(self.nalt)
329	1	alts_out = np.asarray(self.alts).reshape(self.nalt)
330	1	alts_max_out = np.asarray(self.alts_max).reshape(self.nalt)
331
332	1	lats_min_out = np.asarray(self.lats_min).reshape(self.nlat)
333	1	lats_out = np.asarray(self.lats).reshape(self.nlat)
334	1	lats_max_out = np.asarray(self.lats_max).reshape(self.nlat)
335
336	1	ncf = netcdf_file(filename, 'w', **fmtargs)
337
338	1	if self.version is not None:
339	1	ncf.version = self.version
340	1	if self.data_version is not None:
341	1	ncf.L2_data_version = self.data_version
342		#ncf.creation_time = dt.datetime.utcnow().replace(tzinfo=_UTC).strftime("%a %b %d %Y %H:%M:%S %z (%Z)")
343	1	ncf.creation_time = dt.datetime.utcnow().strftime("%a %b %d %Y %H:%M:%S +00:00 (UTC)")
344	1	ncf.author = self.author
345
346		# create netcdf file
347	1	ncf.createDimension('time', None)
348	1	ncf.createDimension('altitude', self.nalt)
349	1	ncf.createDimension('latitude', self.nlat)
350
351	1	forbit = ncf.createVariable('orbit', np.dtype('int32').char, ('time',))
352	1	ftime = ncf.createVariable('time', np.dtype('int32').char, ('time',))
353
354	1	falts_min = ncf.createVariable('alt_min', np.dtype('float64').char, ('altitude',))
355	1	falts = ncf.createVariable('altitude', np.dtype('float64').char, ('altitude',))
356	1	falts_max = ncf.createVariable('alt_max', np.dtype('float64').char, ('altitude',))
357	1	flats_min = ncf.createVariable('lat_min', np.dtype('float64').char, ('latitude',))
358	1	flats = ncf.createVariable('latitude', np.dtype('float64').char, ('latitude',))
359	1	flats_max = ncf.createVariable('lat_max', np.dtype('float64').char, ('latitude',))
360
361	1	falts_min.units = 'km'
362	1	falts_min.positive = 'up'
363	1	falts.units = 'km'
364	1	falts.positive = 'up'
365	1	falts_max.units = 'km'
366	1	falts_max.positive = 'up'
367	1	flats_min.units = 'degrees_north'
368	1	flats.units = 'degrees_north'
369	1	flats_max.units = 'degrees_north'
370
371	1	forbit.units = '1'
372	1	forbit.long_name = 'SCIAMACHY/Envisat orbit number'
373	1	ftime.units = 'days since {0}'.format(self.date0.isoformat(sep=' '))
374	1	ftime.standard_name = 'time'
375
376	1	fdens = ncf.createVariable('density', np.dtype('float64').char, ('time', 'latitude', 'altitude'))
377	1	fdens.units = 'cm^{-3}'
378	1	fdens.standard_name = 'number_concentration_of_nitrogen_monoxide_molecules_in_air'
379	1	fdens_err_meas = ncf.createVariable('error_meas', np.dtype('float64').char, ('time', 'latitude', 'altitude'))
380	1	fdens_err_meas.units = 'cm^{-3}'
381	1	fdens_err_meas.long_name = 'NO number density measurement error'
382	1	fdens_err_tot = ncf.createVariable('error_tot', np.dtype('float64').char, ('time', 'latitude', 'altitude'))
383	1	fdens_err_tot.units = 'cm^{-3}'
384	1	fdens_err_tot.long_name = 'NO number density total error'
385	1	fdens_tot = ncf.createVariable('density_air', np.dtype('float64').char, ('time', 'latitude', 'altitude'))
386	1	fdens_tot.units = 'cm^{-3}'
387	1	fdens_tot.long_name = 'approximate overall number concentration of air molecules (NRLMSIS-00)'
388
389	1	ftime[:] = [self.date]
390	1	forbit[:] = [self.orbit]
391
392	1	falts_min[:] = alts_min_out
393	1	falts[:] = alts_out
394	1	falts_max[:] = alts_max_out
395	1	flats_min[:] = lats_min_out
396	1	flats[:] = lats_out
397	1	flats_max[:] = lats_max_out
398		# reorder by latitude first, then altitude
399	1	fdens[0, :] = self.densities
400		# reorder by latitude first, then altitude
401	1	fdens_err_meas[0, :] = self.dens_err_meas
402	1	fdens_err_tot[0, :] = self.dens_err_tot
403	1	fdens_tot[0, :] = self.dens_tot
404
405		# longitudes if they are available
406	1	if self.nlon > 0:
407	1	lons_out = np.asarray(self.lons).reshape(self.nlon)
408	1	flons = ncf.createVariable('longitude', np.dtype('float64').char, ('time', 'latitude',))
409	1	flons.units = 'degrees_east'
410	1	flons[0, :] = lons_out
411
412	1	if self.apriori is not None:
413	1	fapriori = ncf.createVariable('apriori',
414		np.dtype('float64').char, ('time', 'latitude', 'altitude'))
415	1	fapriori.units = 'cm^{-3}'
416	1	fapriori.long_name = 'apriori NO number density'
417	1	fapriori[0, :] = self.apriori
418
419	1	if self.akdiag is not None:
420	1	fakdiag = ncf.createVariable('akm_diagonal',
421		np.dtype('float64').char, ('time', 'latitude', 'altitude'))
422	1	fakdiag.units = '1'
423	1	fakdiag.long_name = 'averaging kernel matrix diagonal element'
424	1	fakdiag[0, :] = self.akdiag
425	1	if close:
426	1	ncf.close()
427		else:
428	1	return ncf
429
430	1	def read_from_netcdf(self, filename, close=True):
431		"""Read NO densities from netcdf files
432
433		This function has no stream, i.e. file object support.
434
435		Parameters
436		----------
437		filename: str
438		The filename to read the data from.
439		close: bool, optional
440		Whether or not to close the file after reading.
441		Setting to `False` enables reading further data
442		from the same file.
443		Default: True
444
445		Returns
446		-------
447		Nothing if `close` is True. If `close` is False, returns either an
448		`netCDF4.Dataset`,
449		`scipy.io.netcdf.netcdf_file` or
450		`pupynere.netcdf_file` instance depending on availability.
451		"""
452	1	def _try_decode(s):
453	1	if hasattr(s, "decode"):
454		return s.decode()
455	1	return s
456	1	ncf = netcdf_file(filename, 'r')
457
458	1	try:
459	1	self.author = _try_decode(ncf.author)
460		except AttributeError:
461		pass
462	1	try:
463	1	self.version = _try_decode(ncf.version)
464	1	except AttributeError:
465	1	pass
466	1	try:
467	1	self.data_version = _try_decode(ncf.L2_data_version)
468		except AttributeError:
469		pass
470
471	1	self.alts_min = ncf.variables['alt_min'][:].copy()
472	1	self.alts = ncf.variables['altitude'][:].copy()
473	1	self.alts_max = ncf.variables['alt_max'][:].copy()
474	1	self.lats_min = ncf.variables['lat_min'][:].copy()
475	1	self.lats = ncf.variables['latitude'][:].copy()
476	1	self.lats_max = ncf.variables['lat_max'][:].copy()
477
478	1	self.nalt = len(self.alts)
479	1	self.nlat = len(self.lats)
480
481	1	self.date = ncf.variables['time'][:].copy()
482	1	self.orbit = ncf.variables['orbit'][:].copy()
483
484	1	self.densities = ncf.variables['density'][:].copy()
485	1	self.dens_err_meas = ncf.variables['error_meas'][:].copy()
486	1	self.dens_err_tot = ncf.variables['error_tot'][:].copy()
487	1	self.dens_tot = ncf.variables['density_air'][:].copy()
488
489		# longitudes if they are available
490	1	try:
491	1	self.lons = ncf.variables['longitude'][:].copy()
492	1	self.nlon = self.lons.shape[1]
493	1	except KeyError:
494	1	pass
495
496		# apriori
497	1	try:
498	1	self.apriori = ncf.variables['apriori'][:].copy()
499	1	except KeyError:
500	1	pass
501
502		# akm diagonal elements
503	1	try:
504	1	self.akdiag = ncf.variables['akm_diagonal'][:].copy()
505	1	except KeyError:
506	1	pass
507
508	1	if close:
509	1	ncf.close()
510		else:
511		return ncf
512
513	1	def read_from_file(self, filename):
514		"""Wrapper to read NO desnities from files
515
516		Simple wrapper to delegate reading the data from either netcdf
517		or ascii files. Poor man's logic: simply try netcdf first, and
518		if that fails, read as ascii.
519
520		Parameters
521		----------
522		filename: str
523		The filename to read the data from.
524		"""
525	1	try:
526		# try netcdf first
527	1	self.read_from_netcdf(filename)
528	1	except (IOError, OSError, TypeError):
529		# fall back to text file as a last resort
530	1	self.read_from_textfile(filename)
531
532
533	1	def main(*args):
534		argc = len(sys.argv)
535		if argc < 2:
536		print("Not enough arguments, Usage:\n"
537		"{0} [input] output [< input]".format(sys.argv[0]))
538		sys.exit(1)
539		elif argc < 3:
540		try:
541		infile = sys.stdin.buffer # Python 3
542		except AttributeError:
543		infile = sys.stdin
544		outfile = sys.argv[1]
545		else:
546		infile = sys.argv[1]
547		outfile = sys.argv[2]
548		sdl = scia_densities()
549		sdl.read_from_file(infile)
		0 ignored issues – show introduced 2018-07-02 09:55 UTC by Report Bug Copy Issue Report The variable `infile` does not seem to be defined for all execution paths. Loading history...
550		sdl.write_to_netcdf(outfile)
		0 ignored issues – show introduced 2018-07-02 09:55 UTC by Report Bug Copy Issue Report The variable `outfile` does not seem to be defined for all execution paths. Loading history...
551
552
553	1	if __name__ == "__main__":
554		sys.exit(main())
555

st-bender / sciapy

sciapy.level2.density.main() A last analyzed 2023-01-30 12:22 UTC

Complexity

Size

Duplication

Code Coverage

Importance

Duplication Side-by-Side

Filter issues like

sciapy.level2.density.main() A
last analyzed 2023-01-30 12:22 UTC