Updated formatting. (fcd3c6c6) · Commits · CAMS / CSO

py/cso_s5p.py

+55 −44

Original line number	Diff line number	Diff line
		@@ -209,7 +209,6 @@ class S5p_File(object):
		# *

		def GetDim(self, dimname):

		"""
		Return length of dimension.
		"""
		@@ -219,7 +218,6 @@ class S5p_File(object):

		# open:
		with xarray.open_dataset(self.filename, group="PRODUCT") as ds:

		# check ..
		if dimname not in ds.dims:
		logging.error(f"retrieval dimension '{dimname}' not found in: {self.filename}")
		@@ -519,7 +517,9 @@ class S5p_File(object):
		# check units ...
		funits = rcf.get(fkey + ".units")
		if funits != vunits:
		logging.error(f"variable '{varname}' has units '{vunits}' while threshold has units '{funits}'")
		logging.error(
		f"variable '{varname}' has units '{vunits}' while threshold has units '{funits}'"
		)
		raise Exception
		# endif
		# avoid warnings about comparing nan ...
		@@ -541,7 +541,9 @@ class S5p_File(object):
		# check units ...
		funits = rcf.get(fkey + ".units")
		if funits != vunits:
		logging.error(f"variable '{varname}' has units '{vunits}' while threshold has units '{funits}'")
		logging.error(
		f"variable '{varname}' has units '{vunits}' while threshold has units '{funits}'"
		)
		logging.error(f"value range : {numpy.nanmin(values)} - {numpy.nanmax(values)}")
		raise Exception
		# endif
		@@ -565,7 +567,9 @@ class S5p_File(object):
		# check units ...
		funits = rcf.get(fkey + ".units")
		if funits != vunits:
		logging.error(f"variable '{varname}' has units '{vunits}' while threshold has units '{funits}'")
		logging.error(
		f"variable '{varname}' has units '{vunits}' while threshold has units '{funits}'"
		)
		raise Exception
		# endif
		# avoid warnings about comparing nan ...
		@@ -599,11 +603,15 @@ class S5p_File(object):
		# check units ...
		funits = rcf.get(fkey + ".units")
		if funits != vunits:
		logging.error(f"variable '{varname}' has units '{vunits}' while threshold has units '{funits}'")
		logging.error(
		f"variable '{varname}' has units '{vunits}' while threshold has units '{funits}'"
		)
		raise Exception
		# endif
		if funits != vunits2:
		logging.error(f"variable '{varname2}' has units '{vunits2}' while threshold has units '{funits}'")
		logging.error(
		f"variable '{varname2}' has units '{vunits2}' while threshold has units '{funits}'"
		)
		raise Exception
		# endif
		# avoid warnings about comparing nan ...
		@@ -1016,7 +1024,6 @@ class CSO_S5p_File(cso_file.CSO_File):

		# extract 2D array of pixels from 3D track:
		elif oda.dims in [("time", "scanline", "ground_pixel", "level")]:

		# create 2D array with data from selected pixels, copy attributes:
		da = xarray.DataArray(
		numpy.array(
		@@ -1090,7 +1097,6 @@ class CSO_S5p_File(cso_file.CSO_File):

		# extract 3D array of pixels:
		elif oda.dims in [("time", "scanline", "ground_pixel", "level", "level")]:

		# check ...
		if sfile.nlayer is None:
		logging.error("no layer dimension defined yet")
		@@ -1108,7 +1114,9 @@ class CSO_S5p_File(cso_file.CSO_File):
		self.npixel, self.nretr, sfile.nlayer
		),
		dtype=dtype,
		), axes=(0,2,1) ),
		),
		axes=(0, 2, 1),
		),
		dims=vdims,
		coords={"pixel": pixel},
		attrs=oda.attrs,
		@@ -1876,10 +1884,11 @@ class CSO_S5p_File(cso_file.CSO_File):

		# check target dims:
		if vdims == ("pixel", "layeri"):

		# extract per pixel:
		if odat["pmid"].dims in [("time", "scanline", "ground_pixel", "layer"),
		("time", "scanline", "ground_pixel", "level")]:
		if odat["pmid"].dims in [
		("time", "scanline", "ground_pixel", "layer"),
		("time", "scanline", "ground_pixel", "level"),
		]:
		# extract pixels:
		pmid = odat["pmid"].data[itime, iis, iip, :]
		else:
		@@ -1956,10 +1965,11 @@ class CSO_S5p_File(cso_file.CSO_File):

		# check target dims:
		if vdims == ("pixel", "layeri"):

		# extract per pixel:
		if odat["amid"].dims in [("time", "scanline", "ground_pixel", "layer"),
		("time", "scanline", "ground_pixel", "level")]:
		if odat["amid"].dims in [
		("time", "scanline", "ground_pixel", "layer"),
		("time", "scanline", "ground_pixel", "level"),
		]:
		# extract pixels:
		amid = odat["amid"].data[itime, iis, iip, :]
		else:
		@@ -2371,7 +2381,6 @@ class CSO_S5p_File(cso_file.CSO_File):

		# check target dims:
		if vdims == ("pixel", "retr", "retr0"):

		# storage for extract along pixels index:
		pdat = {}
		# extract diagonal:
		@@ -2379,7 +2388,9 @@ class CSO_S5p_File(cso_file.CSO_File):
		if odat[key].dims == ("time", "scanline", "ground_pixel", "level"):
		pdat[key] = odat[key].data[itime, iis, iip, :]
		else:
		logging.error(f"unsupported dimnames '{odat[key].dims}' for special '{special}")
		logging.error(
		f"unsupported dimnames '{odat[key].dims}' for special '{special}"
		)
		raise Exception
		# endif

py/cso_tools.py

+119 −93

Original line number	Diff line number	Diff line
		@@ -33,7 +33,6 @@ Methods


		def linearize_avg_kernel(AK, xa, xa_ratio_max):

		"""
		Linearizes the averaging kernel.

		@@ -83,6 +82,7 @@ def linearize_avg_kernel(AK, xa, xa_ratio_max):

		return AK_lin


		# enddef linearize_avg_kernel


		@@ -193,7 +193,7 @@ def compute_ellipse_footprint_array(geodeticLat, longitude, range_sat, azimuth,
		losPointX3 = rotate_vector_array(rotationVector, geodeticLatRad, losPointX2)
		rotationVector = np.array([0.0e0, 0.0e0, 1.0e0])
		useLonRad = lonRad.copy()
		selb, = np.where(lonRad < 0.0)
		(selb,) = np.where(lonRad < 0.0)
		if len(selb) >= 0:
		useLonRad[selb] = lonRad[selb] + 2.0 * np.pi
		# else:
		@@ -236,7 +236,9 @@ def compute_ellipse_footprint_array(geodeticLat, longitude, range_sat, azimuth,
		# Find the geodetic latitude and longitude ( geolocation)
		# for the given current fov Vector and satellite position.
		#
		geolocationLat, geolocationLon = compute_geolocation_array(satellitePosition, curFovVector, earthRadius, flatFact)
		geolocationLat, geolocationLon = compute_geolocation_array(
		satellitePosition, curFovVector, earthRadius, flatFact
		)

		# format long eng
		# satellitePosition
		@@ -249,19 +251,17 @@ def compute_ellipse_footprint_array(geodeticLat, longitude, range_sat, azimuth,
		ellLat[:, iAngleStep] = geolocationLat
		ellLon[:, iAngleStep] = geolocationLon


		# end


		return ellLat, ellLon


		# endef compute_ellipse_footprint_array


		# *



		def compute_geolocation(satellitePosition, lineOfSight, earthRadius, flatFact):
		#
		# Description:
		@@ -333,32 +333,41 @@ def compute_geolocation(satellitePosition, lineOfSight, earthRadius, flatFact):
		# a quadratic equation where A lambda**2 + B lambda + C = 0, here x is the slant range.
		#
		polarRadius = earthRadius * (1.0 - flatFact)
		termA = ((lineOfSight[0] / earthRadius) 2) + ((lineOfSight[1] / earthRadius) 2) + (
		(lineOfSight[2] / polarRadius) ** 2)

		termB = (satellitePosition[0] * lineOfSight[0] / (earthRadius ** 2)) + (
		satellitePosition[1] * lineOfSight[1] / (earthRadius ** 2)) + (
		satellitePosition[2] * lineOfSight[2] / (polarRadius ** 2))
		termA = (
		((lineOfSight[0] / earthRadius) ** 2)
		+ ((lineOfSight[1] / earthRadius) ** 2)
		+ ((lineOfSight[2] / polarRadius) ** 2)
		)

		termB = (
		(satellitePosition[0] * lineOfSight[0] / (earthRadius**2))
		+ (satellitePosition[1] * lineOfSight[1] / (earthRadius**2))
		+ (satellitePosition[2] * lineOfSight[2] / (polarRadius**2))
		)
		termB = termB * 2.0

		termC = ((satellitePosition[0] / earthRadius) 2) + ((satellitePosition[1] / earthRadius) 2) + (
		(satellitePosition[2] / polarRadius) ** 2) - 1.0
		termC = (
		((satellitePosition[0] / earthRadius) ** 2)
		+ ((satellitePosition[1] / earthRadius) ** 2)
		+ ((satellitePosition[2] / polarRadius) ** 2)
		- 1.0
		)

		radical = termB*2 - (4.0 termA * termC)

		if (radical < 0.0):
		if radical < 0.0:
		# The line of sight does not intercept the Earth ellipsoid.
		#
		ellLat = -999.0
		ellLon = -999.0
		# end
		if (radical == 0.0):
		if radical == 0.0:
		# The line of sight does not intercept the Earth ellipsoid tangentially.
		#
		slantRange = -termB / (2.0 * termA)
		geolocationPoint = satellitePosition + slantRange * lineOfSight
		# end
		if (radical > 0.0):
		if radical > 0.0:
		# The line of sight intercepts the Earth ellipsoid at 2 point, the solution
		# is the shorter slant range.
		#
		@@ -375,13 +384,13 @@ def compute_geolocation(satellitePosition, lineOfSight, earthRadius, flatFact):

		return geoLat, geoLon


		# endef compute_geolocation


		# *



		def compute_geolocation_array(satellitePosition, lineOfSight, earthRadius, flatFact):
		#
		# Description:
		@@ -453,30 +462,41 @@ def compute_geolocation_array(satellitePosition, lineOfSight, earthRadius, flatF
		# a quadratic equation where A lambda**2 + B lambda + C = 0, here x is the slant range.
		#
		polarRadius = earthRadius * (1.0 - flatFact)
		termA = ((lineOfSight[:,0] / earthRadius) 2) + ((lineOfSight[:,1] / earthRadius) 2) + (
		(lineOfSight[:,2] / polarRadius) ** 2)

		termB = (satellitePosition[:,0] * lineOfSight[:,0] / (earthRadius ** 2)) + (
		satellitePosition[:,1] * lineOfSight[:,1] / (earthRadius ** 2)) + (
		satellitePosition[:,2] * lineOfSight[:,2] / (polarRadius ** 2))
		termA = (
		((lineOfSight[:, 0] / earthRadius) ** 2)
		+ ((lineOfSight[:, 1] / earthRadius) ** 2)
		+ ((lineOfSight[:, 2] / polarRadius) ** 2)
		)

		termB = (
		(satellitePosition[:, 0] * lineOfSight[:, 0] / (earthRadius**2))
		+ (satellitePosition[:, 1] * lineOfSight[:, 1] / (earthRadius**2))
		+ (satellitePosition[:, 2] * lineOfSight[:, 2] / (polarRadius**2))
		)
		termB = termB * 2.0

		termC = ((satellitePosition[:,0] / earthRadius) 2) + ((satellitePosition[:,1] / earthRadius) 2) + (
		(satellitePosition[:,2] / polarRadius) ** 2) - 1.0
		termC = (
		((satellitePosition[:, 0] / earthRadius) ** 2)
		+ ((satellitePosition[:, 1] / earthRadius) ** 2)
		+ ((satellitePosition[:, 2] / polarRadius) ** 2)
		- 1.0
		)

		radical = termB*2 - (4.0 termA * termC)

		selb, = np.where(radical<0.0)
		sele, = np.where(radical==0.0)
		sela, = np.where(radical>0.0)
		if len(sele)>0.:
		(selb,) = np.where(radical < 0.0)
		(sele,) = np.where(radical == 0.0)
		(sela,) = np.where(radical > 0.0)
		if len(sele) > 0.0:
		slantRange = -termB[sele] / (2.0 * termA[sele])
		geolocationPoint[sele] = satellitePosition[sele] + slantRange * lineOfSight[sele]
		if len(sela) > 0:
		slantRange1 = (-termB[sela] - np.sqrt(radical[sela])) / (2.0 * termA[sela])
		slantRange2 = (-termB[sela] + np.sqrt(radical[sela])) / (2.0 * termA[sela])
		slantRange = np.min([slantRange1, slantRange2], 0)
		geolocationPoint[sela,:] = satellitePosition[sela,:] + slantRange[:,None] * lineOfSight[sela,:]
		geolocationPoint[sela, :] = (
		satellitePosition[sela, :] + slantRange[:, None] * lineOfSight[sela, :]
		)
		# if (radical < 0.0):
		# # The line of sight does not intercept the Earth ellipsoid.
		# #
		@@ -506,6 +526,7 @@ def compute_geolocation_array(satellitePosition, lineOfSight, earthRadius, flatF

		return geoLat, geoLon


		# enddef compute_geolocation_array


		@@ -556,13 +577,14 @@ def geocentric_range(geocentricLatitude, equatorialRadius, flatteningFactor):
		# modules:
		import numpy as np

		deg2rad = np.pi / 180.0;
		deg2rad = np.pi / 180.0
		polarRadius = (1.0 - flatteningFactor) * equatorialRadius
		lat = geocentricLatitude * deg2rad
		range = np.sqrt((equatorialRadius * np.cos(lat)) ** 2 + (polarRadius * np.sin(lat)) ** 2)

		return range


		# enddef geocentric_range


		@@ -622,12 +644,10 @@ def latitude_transform_array(latRef, flatteningFactor, option):
		latIn = latRef * deg2rad
		ff = 1.0 - flatteningFactor

		if (option == 1):

		if option == 1:
		latOut = np.arctan(((ff) ** 2) * np.tan(latIn))
		latitudeOutput = latOut / deg2rad
		else:

		latOut = np.arctan((1.0 / (ff*2)) np.tan(latIn))
		latitudeOutput = latOut / deg2rad

		@@ -635,6 +655,7 @@ def latitude_transform_array(latRef, flatteningFactor, option):

		return latitudeOutput


		# enddef latitude_transform_array


		@@ -698,10 +719,13 @@ def rotate_vector_array(rotationAxis, angle, oldVector):

		firstTerm = oldVector * cosAngle[:, None]
		secondTerm = np.cross(oldVector, rotationAxis) * sinAngle[:, None]
		thirdTerm = (np.dot(oldVector, rotationAxis) * rotationAxis[:,None]).T * (1.0 - cosAngle)[:,None]
		thirdTerm = (np.dot(oldVector, rotationAxis) * rotationAxis[:, None]).T * (1.0 - cosAngle)[
		:, None
		]
		newVector = firstTerm - secondTerm + thirdTerm
		return newVector


		# enddef rotate_vector_array


		@@ -765,12 +789,14 @@ def rotate_vector_array_adjusted(rotationAxis, angle, oldVector):

		firstTerm = oldVector * cosAngle[:, None]
		secondTerm = np.cross(oldVector, rotationAxis) * sinAngle[:, None]
		thirdTerm = (np.sum(oldVector* rotationAxis,1)[:,None] * rotationAxis) * (1.0 - cosAngle)[:,None]
		thirdTerm = (np.sum(oldVector * rotationAxis, 1)[:, None] * rotationAxis) * (1.0 - cosAngle)[
		:, None
		]
		newVector = firstTerm - secondTerm + thirdTerm
		return newVector

		# enddef rotate_vector_array_adjusted

		# enddef rotate_vector_array_adjusted


		########################################################################

py/cso_dataspace.py

+1 −1

File changed.

Contains only whitespace changes.

py/cso_plot.py

+1 −1

File changed.

Contains only whitespace changes.