sciapy.level2.density.scia_densities.__init__() - Code Metrics - Inspection of "ci: Fix builds on travis" - st-bender/sciapy - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( 52d259...18360e )

by Stefan

created 2020-06-15 20:27 UTC

sciapy.level2.density.scia_densities.init() A

↳ Parent: sciapy.level2.density

Complexity

Conditions

Size

Total Lines	19
Code Lines	19

Duplication

Lines	19
Ratio	100 %

Code Coverage

Tests	19
CRAP Score	1

Importance

Changes

Metric	Value
cc	1
eloc	19
nop	5
dl	19
loc	19
ccs	19
cts	19
cp	1
crap	1
rs	9.45
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	View Code Duplication	class scia_densities(object):
			0 ignored issues – show Duplication introduced 2020-03-03 10:24 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
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	View Code Duplication	def main(*args):
			0 ignored issues – show Duplication introduced 2018-07-02 09:55 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
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

Push — master ( 52d259...18360e )

sciapy.level2.density.scia_densities.__init__() A

Complexity

Size

Duplication

Code Coverage

Importance

Duplication Side-by-Side

Filter issues like

sciapy.level2.density.scia_densities.init() A