Support orbits that overlap with dateline. (33a3b05a) · Commits · CAMS / CSO

py/cso_gridded.py

+18 −0

Original line number	Diff line number	Diff line
		@@ -23,6 +23,10 @@
		# 2023-08, Arjo Segers
		# Define lon/lat axes as dtype f4.
		#
		# 2024-01, Arjo Segers
		# Support pixels that overlap with dateline and have corners with
		# longitudes outside [-180,180].
		#


		########################################################################
		@@ -678,6 +682,20 @@ class CSO_GriddedAverage(utopya.UtopyaRc):
		ii2 = df2.index.get_level_values("ii").values
		jj2 = df2.index.get_level_values("jj").values

		# in orbits over the date-line the corners of the footprints
		# might have longitudes outside [-180,+180] to maintain a convex region ...
		# check if grid spans the globe in longitudes (east == west):
		if abs(east - west) % 360.0 < 1.0e-3:
		# reset values left from domain:
		(kk,) = numpy.where(ii2 < 0)
		if len(kk) > 0:
		ii2[kk] += nlon
		# also for values right from domain:
		(kk,) = numpy.where(ii2 >= nlon)
		if len(kk) > 0:
		ii2[kk] -= nlon
		# endif # global grid

		# target array for total pixel weight per grid cell:
		pat_ww = numpy.zeros((nlat, nlon), float)
		# within grid:

py/cso_s5p.py

+124 −5

Original line number	Diff line number	Diff line
		@@ -35,6 +35,11 @@
		# 2024-01, Arjo Segers
		# Switch between DataSpace and PAL downloader based on download link.
		#
		# 2024-01, Arjo Segers
		# For pixels over dateline, reset longitude bounds to values
		# outside [-180,+180] to define the footprint as a convex shape
		# that includes the pixel center.
		#

		########################################################################
		###
		@@ -638,10 +643,19 @@ class CSO_S5p_File(cso_file.CSO_File):
		For these variables a key '``special``' is used which will enable the
		correct conversion. The following specials are currently implemented:

		* ``track_longitude`` : longiudes at centers of original 2D track; requires a ``.from`` setting
		* ``track_latitude`` : latiudes at centers of original 2D track; requires a ``.from`` setting
		* ``track_longitude_bounds`` : longiude bounds at centers of original 2D track; requires a ``.from`` setting
		* ``track_latitude_bounds`` : latiude bounds at centers of original 2D track; requires a ``.from`` setting
		* ``longitude_bounds`` : longitude bounds per pixel;
		if pixel covers date line, corner values are reset to values outside [-180,+180]
		to ensure that the footprints remains convex with the center inside::

		<rcbase>.output.var.longitude_bounds.from : PRODUCT/SUPPORT_DATA/GEOLOCATIONS/longitude_bounds
		<rcbase>.output.var.longitude_bounds.special : longitude_bounds

		* ``track_longitude`` : longitudes at centers of original 2D track; requires a ``.from`` setting
		* ``track_latitude`` : latitudes at centers of original 2D track; requires a ``.from`` setting
		* ``track_longitude_bounds`` : longitude bounds of original 2D track;
		if pixel covers date line, corner values are reset to values outside [-180,+180]
		to ensure that the footprints remains convex with the center inside;
		* ``track_latitude_bounds`` : latitude bounds of original 2D track; requires a ``.from`` setting
		* ``ground_pixel`` : index of ground pixel in original 2D track; requires a ``.from`` setting
		* ``sum`` : create a variable as the sum over over layers; requires a ``.from`` setting
		* ``square`` : create a variable as the square of the input; requires a ``.from`` setting
		@@ -1220,6 +1234,70 @@ class CSO_S5p_File(cso_file.CSO_File):
		da.attrs["long_name"] = "time"
		da.encoding["calendar"] = "standard"

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# pixel corners: for longitudes, need to ensure that also near dateline
		# the corners form a convex region around center
		elif special in ["longitude_bounds"]:
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# name of original field:
		opath = rcf.get(vkey + ".from")
		# get data array:
		oda = sfile.GetData(opath, units=vunits)

		# info ..
		logging.info(" wrap longitudes around the world ...")
		# reset longitudes to be wrapped around the world in one contineous series:
		oda__data = self.ResetLongitudeBounds(oda.data[itime, :])
		# for list of pixels the values could be reset to [-180,180] again if not crossing the date line ..
		# ~ left from dateline?
		while True:
		# select pixels with all corners left from dateline:
		inds = numpy.where(numpy.all(oda__data < -180.0, axis=2))
		# no pixels found? leave loop:
		if len(inds[0]) == 0:
		break
		# reset:
		oda__data[inds] += 360.0
		# endwhile
		# ~ right from dateline?
		while True:
		# select pixels with all corners right from dateline:
		inds = numpy.where(numpy.all(oda__data > 180.0, axis=2))
		# no pixels found? leave loop:
		if len(inds[0]) == 0:
		break
		# reset:
		oda__data[inds] -= 360.0
		# endwhile

		# ~ corners per pixel:
		if vdims == ("pixel", "corner"):

		# bound arrays for longitude etc:
		if oda.dims in [("time", "scanline", "ground_pixel", "corner")]:

		# create 1D array with data from selected pixels, copy attributes:
		da = xarray.DataArray(
		numpy.array(oda__data[iis, iip, :], dtype=dtype),
		dims=vdims,
		coords={"pixel": pixel},
		attrs=oda.attrs,
		)

		else:
		logging.error(
		f"could not convert from dims '{oda.dims}'"
		+ " to target dims '{vdims}'"
		)
		raise Exception
		# endif
		else:
		logging.error(f"unsupported target dims: {vdims}")
		logging.error(f" target variable : {varname}")
		raise Exception
		# endif

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# full track, centers:
		elif special in ["track_longitude", "track_latitude"]:
		@@ -1259,7 +1337,48 @@ class CSO_S5p_File(cso_file.CSO_File):

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# full track, corners:
		elif special in ["track_longitude_bounds", "track_latitude_bounds"]:
		elif special in ["track_longitude_bounds"]:
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# name of original field:
		opath = rcf.get(vkey + ".from")
		# get data array:
		oda = sfile.GetData(opath, units=vunits)

		# info ..
		logging.info(" wrap longitudes around the world ...")
		# reset longitudes to be wrapped around the world in one contineous series:
		oda__data = self.ResetLongitudeBounds(oda.data[itime, is0:is1, :])

		# check target dims:
		if vdims == ("track_scan", "track_pixel", "corner"):

		# check shape:
		if oda.dims in [("time", "scanline", "ground_pixel", "corner")]:

		# create 2D array (+corners) with data from selected scan lines:
		da = xarray.DataArray(
		oda__data,
		dims=vdims,
		coords={"track_scan": track_scan, "track_pixel": track_pixel,},
		attrs=oda.attrs,
		)
		else:
		logging.error(
		'could not convert from dims "%s" to target dims "%s"'
		% (str(oda.dims), str(vdims))
		)
		raise Exception
		# endif
		else:
		logging.error("unsupported target dims: %s" % str(vdims))
		logging.error(" target variable : %s" % varname)
		raise Exception
		# endif

		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# full track, corners:
		elif special in ["track_latitude_bounds"]:
		# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		# name of original field:

py/cso_superobs.py

+18 −0

Original line number	Diff line number	Diff line
		@@ -11,6 +11,10 @@
		# 2023-01, Arjo Segers
		# Fixed averaging over time values.
		#
		# 2024-01, Arjo Segers
		# Support pixels that overlap with dateline.
		#


		########################################################################
		###
		@@ -508,6 +512,20 @@ class CSO_SuperObs(utopya.UtopyaRc):
		# idem for time:
		df2__time = df__time.groupby(["jj", "ii"]).mean()

		# in orbits over the date-line the corners of the footprints
		# might have longitudes outside [-180,+180] to maintain a convex region ...
		# check if grid spans the globe in longitudes (east == west):
		if abs(east - west) % 360.0 < 1.0e-3:
		# reset values left from domain:
		(kk,) = numpy.where(ii2 < 0)
		if len(kk) > 0:
		ii2[kk] += nlon
		# also for values right from domain:
		(kk,) = numpy.where(ii2 >= nlon)
		if len(kk) > 0:
		ii2[kk] -= nlon
		# endif # global grid

		# mask for super pixels with enough contributions(by number or area)
		# of original pixels:
		if wtype == "number":