Plot pixels als polygons in case track information is not present. (436d7701) · Commits · CAMS / CSO

src/cso/cso_catalogue.py

+131 −47

Original line number	Diff line number	Diff line
		@@ -35,6 +35,11 @@
		# Extended unit conversions.
		# Fixed reading of settings for map styles.
		#
		# 2025-09, Arjo Segers
		# Plot pixels als polygons in case track information is not present.
		# Extended unit conversions.
		# Fixed string formatting.
		#


		########################################################################
		@@ -156,6 +161,30 @@ class CSO_CatalogueBase(utopya.UtopyaRc):
		# apply factor:
		values = values * 1e6

		# vertical column density:
		elif conversion == "molec/cm2 -> umol/m2":
		# apply factor:
		# molec/cm2 / (molec/mol) umol/mol cm2/m2
		values = values / cso_constants.Avog * 1e6 * 1e4

		# vertical column density:
		elif conversion == "molec/cm2 -> 1e15 mlc/cm2":
		# apply factor:
		# molec/cm2
		values = values / 1e15

		# vertical column density:
		elif conversion == "(molec/cm2)**2 -> 1e15 mlc/cm2":
		# apply factor:
		# molec/cm2
		values = numpy.sqrt(values) / 1e15

		# vertical column density, squared
		elif conversion == "(molec/cm2)**2 -> umol/m2":
		# apply factor:
		# molec/cm2 / (molec/mol) umol/mol cm2/m2
		values = numpy.sqrt(values) / cso_constants.Avog * 1e6 * 1e4

		# vertical column density:
		elif conversion == "mol m-2 -> mmol/m2":
		# apply factor:
		@@ -383,7 +412,11 @@ class CSO_Catalogue(CSO_CatalogueBase):
		filename = os.path.join(inputdir, fname)

		# orbit key:
		if "orbit" in xlst:
		orbit = xlst.GetValue(fname, "orbit")
		else:
		orbit = ""
		#endif

		## testing ...
		# if '1200' not in filename : continue
		@@ -444,23 +477,32 @@ class CSO_Catalogue(CSO_CatalogueBase):
		# settings for this variable:
		vkey = "var.%s" % varname
		# plot type:
		ptype = self.GetSetting("var.{varname}.type", default="map")
		ptype = self.GetSetting(f"var.{varname}.type", default="map")
		# input variable:
		vname = self.GetSetting("var.{varname}.source", default=varname)
		vname = self.GetSetting(f"var.{varname}.source", default=varname)
		# target units:
		vunits = self.GetSetting("var.{varname}.units", default="None")
		vunits = self.GetSetting(f"var.{varname}.units", default="None")
		# long name used in labels:
		long_name = self.GetSetting("var.{varname}.long_name", default=varname)
		long_name = self.GetSetting(f"var.{varname}.long_name", default=varname)

		# track defined in files?
		on_track = "track_longitude" in orb.ds

		# switch:
		if ptype == "map":
		# style:
		vmin = eval(self.GetSetting("var.{varname}.vmin", default="None"))
		vmax = eval(self.GetSetting("var.{varname}.vmax", default="None"))
		colors = eval(self.GetSetting("var.{varname}.colors", default="None"))
		vmin = eval(self.GetSetting(f"var.{varname}.vmin", default="None"))
		vmax = eval(self.GetSetting(f"var.{varname}.vmax", default="None"))
		colors = eval(self.GetSetting(f"var.{varname}.colors", default="None"))

		# extract corner grids and values:
		xx, yy, values, attrs = orb.GetTrack(vname)
		if on_track:
		xx, yy, values, attrs = orb.GetTrack(sstate.ds[vname])
		else:
		xx, attrs = orb.GetData( "longitude_bounds" )
		yy, attrs = orb.GetData( "latitude_bounds" )
		values, attrs = orb.GetData( vname )
		#endif

		# adhoc, some chocho files have nans in lon/lat bound arrays,
		# these are still present because the scan line has some valid values ...
		@@ -512,6 +554,7 @@ class CSO_Catalogue(CSO_CatalogueBase):
		bmp.update(bmp_kwargs)

		# create map figure:
		if on_track:
		fig = cso_plot.QuickPat(
		values,
		x=xx,
		@@ -525,6 +568,21 @@ class CSO_Catalogue(CSO_CatalogueBase):
		title=title,
		figsize=figsize,
		)
		else:
		fig = cso_plot.QuickPoly(
		values,
		x=xx,
		y=yy,
		vmin=vmin,
		vmax=vmax,
		cmap=dict(colors=colors, color_bad=color_nan),
		bmp=bmp,
		domain=domain,
		cbar=dict(label=label),
		title=title,
		figsize=figsize,
		)
		#endif
		# save:
		fig.Export(figfile)

		@@ -812,25 +870,35 @@ class CSO_SimCatalogue(CSO_CatalogueBase):
		else:
		vfile, vname = "data", vsource
		# endif
		# track defined in files?
		on_track = "track_longitude" in sdata.ds
		# extract corner grids and values:
		if vfile == "state":
		# check ...
		if vname not in sstate.ds.keys():
		logging.error(
		"variable `%s` not found in state from file: %s"
		% (vname, statefile)
		)
		logging.error(f"variable '{vname}' not found in state from file: {statefile}")
		raise Exception
		# endif
		# extract corner grids and values:
		if on_track:
		xx, yy, values, attrs = sdata.GetTrack(sstate.ds[vname])
		else:
		xx, attrs = sdata.GetData( "longitude_bounds" )
		yy, attrs = sdata.GetData( "latitude_bounds" )
		values, attrs = sstate.GetData( vname )
		#endif
		elif vfile == "data":
		# extract corner grids and values from data file:
		if on_track:
		xx, yy, values, attrs = sdata.GetTrack(vname)
		else:
		xx, attrs = sdata.GetData( "longitude_bounds" )
		yy, attrs = sdata.GetData( "latitude_bounds" )
		values, attrs = sdata.GetData( vname )
		#endtry
		else:
		logging.error(
		'unsupported type keyword "%s" in source description "%s" for variable "%s"'
		% (vfile, vsource, varname)
		f"unsupported type keyword '{vfile}' in source description '{vsource}' for variable '{varname}'"
		)
		raise Exception
		# endif
		@@ -849,9 +917,7 @@ class CSO_SimCatalogue(CSO_CatalogueBase):
		if (len(shp) == 3) and (shp[2] == 1):
		values = values[:, :, 0]
		else:
		logging.error(
		'could not squeeze "%s" %s to 2D track' % (vname, str(shp))
		)
		logging.error(f"could not squeeze '{vname}' {shp} to 2D track")
		raise Exception
		# endif
		# endif
		@@ -866,7 +932,10 @@ class CSO_SimCatalogue(CSO_CatalogueBase):
		bmp = dict(countries=True)
		bmp.update(bmp_kwargs)

		# info ...
		logging.info(indent + " plot pattern ...")
		# create map figure:
		if on_track:
		fig = cso_plot.QuickPat(
		values,
		x=xx,
		@@ -880,6 +949,21 @@ class CSO_SimCatalogue(CSO_CatalogueBase):
		title=title,
		figsize=figsize,
		)
		else:
		fig = cso_plot.QuickPoly(
		values,
		x=xx,
		y=yy,
		vmin=vmin,
		vmax=vmax,
		cmap=dict(colors=colors, color_bad=color_nan),
		bmp=bmp,
		domain=domain,
		cbar=dict(label=label),
		title=title,
		figsize=figsize,
		)
		#endif
		# save:
		fig.Export(figfile)

		@@ -902,7 +986,7 @@ class CSO_SimCatalogue(CSO_CatalogueBase):
		# break

		else:
		logging.info(indent + " no state file: %s" % statefile)
		logging.info(f"{indent} no state file: {statefile}")
		# endif # state file present

		# next:

src/cso/cso_plot.py

+159 −69

Original line number	Diff line number	Diff line
		@@ -38,6 +38,9 @@
		# Improved supported for defined number of color entries.
		# Support masked arrays for selecting no-data values.
		#
		# 2025-09, Arjo Segers
		# Implemented 'QuickPoly' method to plot list of polygons.
		#

		#
		# NOTES
		@@ -1195,6 +1198,7 @@ class ColorbarFigure(Figure):
		# *

		def Polygons(self, xx, yy, values, vmin=None, vmax=None, vbounds=None, **kwargs):

		"""
		Add colored polygons.

		@@ -1211,81 +1215,36 @@ class ColorbarFigure(Figure):

		"""

		# obtain mappable for colorbar,
		# easiest was just to mark the centers with a tiny dot;
		# this will define the "cmap" and "cbar" attributes:
		p = self.Scatter(
		xx.mean(axis=1),
		yy.mean(axis=1),
		c=values,
		vmin=vmin,
		vmax=vmax,
		vbounds=vbounds,
		marker=".",
		markersize=0,
		linewidths=0,
		)

		# map values to colors:
		colors = self.cmap(self.cnorm(values))
		# loop:
		for ip in range(len(values)):
		# add patch:
		self.ax.fill(xx[ip, :], yy[ip, :], color=colors[ip])
		# endfor

		# enddef Polygons

		# *

		def PolygonsMesh(self, xx, yy, values, vmin=None, vmax=None, vbounds=None, **kwargs):
		"""
		Add colored polygons.

		Arguments:

		* ``xx``, ``yy`` : 2D corner fields, shapes (ny+1,nx+1)

		Optional arguments:

		* ``vmin``,``vmax`` : data range for coloring, min and max values by default,
		or the settings for PColor if that was called before
		* ``vbounds`` : value boundaries (alternative for vmin,vmax)
		# modules:
		import numpy
		import matplotlib.pyplot as plt
		import matplotlib.cm as cm
		import matplotlib.collections

		"""
		# collect vertices into array with shape (npoly,ncorner,2)
		verts = numpy.array( [xx.T,yy.T] ).T

		# to obtain mappable for colorbar,
		# easiest was just to mark the centers with a tiny dot;
		# this will define the "cmap" and "cbar" attributes:
		p = self.Scatter(
		(xx[:-1, :-1] + xx[1:, :-1] + xx[:-1, 1:] + xx[1:, 1:]) / 4.0,
		(yy[:-1, :-1] + yy[1:, :-1] + yy[:-1, 1:] + yy[1:, 1:]) / 4.0,
		c=values,
		vmin=vmin,
		vmax=vmax,
		vbounds=vbounds,
		marker=".",
		linewidths=0,
		)
		# set mapping from data to color if not done yet:
		if self.cnorm is None:
		# set norm function:
		self.Set_ColorNorm(values, vmin=vmin, vmax=vmax, vbounds=vbounds)
		# endif

		# map values to colors:
		colors = self.cmap(self.cnorm(values))

		# shape:
		ny, nx = values.shape
		# loop:
		for j in range(ny):
		for i in range(nx):
		# add patch:
		self.ax.fill(
		[xx[j, i], xx[j, i + 1], xx[j + 1, i + 1], xx[j + 1, i]],
		[yy[j, i], yy[j, i + 1], yy[j + 1, i + 1], yy[j + 1, i]],
		color=colors[j, i],
		)
		# endfor # i
		# endfor # j
		# define polygons:
		collection = matplotlib.collections.PolyCollection( verts,
		facecolors=colors, edgecolors="face",
		cmap=self.cmap, norm=self.cnorm )

		# enddef PolygonsMesh
		# add to figure:
		p = self.ax.add_collection( collection )

		# add colorbar if not done yet:
		self.Add_ColorBar(p)

		# enddef Polygons

		# *

		@@ -1734,6 +1693,137 @@ def QuickPat(
		# *


		def QuickPoly(
		cc,
		x=None,
		y=None,
		domain=None,
		title=None,
		vmin=None,
		vmax=None,
		dmax=None,
		vbounds=None,
		bmp=None,
		cbar=dict(),
		border=None,
		fig=None,
		**kwargs,
		):
		"""

		Plot series of polygons.

		Arguments passed to :py:meth:`ColorbarFigure.Polygons` method:

		* ``x``, ``y`` : polygon coordinates with shape ``(npoly,ncorner)``;
		* ``domain=(left,right,bottom,top)`` : coordinate bounds; if ``domain`` is ``None`` or
		some elements are ``None``, values are based on the ``x`` or ``y`` coordinates if defined.
		* ``cc`` : data values, shape (ny,nx)
		* ``vmin``,``vmax`` : data range for coloring
		* ``vbounds`` : value bounds for colorbar (defines interval with same color)

		Other arguments are passed to 'ColorbarFigure' :

		* ``figsize = (width,height)``
		* ``domain=(left,right,bottom,top)`` : axes extent
		* ``title`` : text above axes
		* ``cmap = dict( 'colors'=['blue','white','red'], ncolor=9 )``
		* ``cbar = dict( 'ntick'=9, extend='both', label='entity [units]' )``,
		or ``False`` to skip colorbar because it is plotted elsewhere.
		* ``bmp`` : map projection, see :py:meth:`.Figure.AddAxes`::

		bmp=dict( countries=True, raster=[30.0,30.0] )

		Return value:

		* ``fig`` : Figure instance

		If an already existing Figure instance is passed as the ``fig`` argument,
		the pattern will be added to the plot.

		"""

		# modules:
		import numpy

		# set domain if not present yet:
		if domain is None:
		domain = [x.min(), x.max(), y.min(), y.max()]
		else:
		if domain[0] is None:
		domain[0] = xx.min()
		if domain[1] is None:
		domain[1] = xx.max()
		if domain[2] is None:
		domain[2] = yy.min()
		if domain[3] is None:
		domain[3] = yy.max()
		# endif

		# replace or add domain in bmp argument:
		if type(bmp) == dict:
		bmp["domain"] = domain

		# symetric range?
		if dmax is not None:
		# range:
		vmin = -dmax
		vmax = dmax
		# colors:
		if "cmap" in kwargs.keys():
		if "colors" not in kwargs["cmap"].keys():
		kwargs["cmap"]["colors"] = "pwb"
		# endif
		else:
		kwargs["cmap"] = {"colors": "pwb"}
		# endif
		# endif

		# value range:
		if hasattr(cc, "values"):
		cmin = numpy.nanmin(cc.values)
		cmax = numpy.nanmax(cc.values)
		else:
		cmin = numpy.nanmin(cc)
		cmax = numpy.nanmax(cc)
		# endif

		# auto extend?
		if "extend" not in cbar.keys():
		# extend below minimum?
		extend_min = (vmin is not None) and (cmin < vmin)
		# extend above maximum?
		extend_max = (vmax is not None) and (cmax > vmax)
		# set keyword:
		if extend_min and (not extend_max):
		cbar["extend"] = "min"
		elif extend_min and extend_max:
		cbar["extend"] = "both"
		elif (not extend_min) and extend_max:
		cbar["extend"] = "max"
		else:
		cbar["extend"] = "neither"
		# endif
		# endif

		# setup figure:
		if fig is None:
		fig = ColorbarFigure(cbar=cbar, bmp=bmp, domain=domain, title=title, **kwargs)
		# endif

		# add polygons:
		fig.Polygons(x, y, cc, vmin=vmin, vmax=vmax, vbounds=vbounds)

		# ok
		return fig


		# enddef QuickPoly


		# *


		def QuickDens(x, y, bins=10, normed=False, **kwargs):
		"""
		Density plot.