Support computation of grid cell areas for 2D arrays. (928beb41) · Commits · CAMS / CSO

py/cso_mapping.py

+59 −24

Original line number	Diff line number	Diff line
		#
		# Changes
		#
		# 2022-09, Arjo Segers, Met-Norway
		# Support computation of grid cell areas for 2D arrays.
		#

		########################################################################
		###
		### help
		###
		########################################################################


		"""
		**********************
		@@ -186,7 +199,7 @@ def LonLatRectanglesArea( xx, yy ) :
		"""
		Return areas in m2 of rectangles defined by:

		* ``xx``, ``yy`` : corner locations in degrees east/north as ``(n,4)`` arrays
		* ``xx``, ``yy`` : corner locations in degrees east/north as ``(n,4)`` or ``(m,n,4)`` arrays

		Approximate area of rectangle in ``(lon,lat)`` in degrees using 2 triangles::

		@@ -201,12 +214,33 @@ def LonLatRectanglesArea( xx, yy ) :
		# modules:
		import numpy

		# switch:
		if (len(xx.shape) == 2) and (len(yy.shape) == 2) :
		# area's of first triangles:
		A1 = LonLatTrianglesArea( xx[:,0:3], yy[:,0:3] )
		# area's of second triangles:
		A2 = LonLatTrianglesArea( numpy.hstack((xx[:,2:4],xx[:,0:1])), numpy.hstack((yy[:,2:4],yy[:,0:1])) )
		# combine:
		return A1 + A2
		A = A1 + A2

		elif (len(xx.shape) == 3) and (len(yy.shape) == 3) :
		# shape:
		m,n,nc = xx.shape
		# init output:
		A = numpy.zeros((m,n),float)
		# loop over rows:
		for j in range(m) :
		# recursive call:
		A[j,:] = LonLatRectanglesArea( xx[j,:], yy[j,:] )
		#endfor

		else :
		logging.error( 'unsupported input shapes xx(%s) and yy(%s)' % (str(xx.shape),str(yy.shape)) )
		raise Exception
		#endif

		# ok:
		return A

		#enddef LonLatRectanglesArea

		@@ -222,7 +256,7 @@ def LonLatPolygonCentroids( xx, yy, maxlevel=5, _level=0, indent='' ) :

		"""
		Let :math:`N` polygons be defined with :math:`n_c` corners in the arguments:
		This method return arrays with :math:`P` locations of points distributed within the polygons,
		This method returns arrays with :math:`(P,N)` locations of points distributed within the polygons,
		as well as the fraction of the area corresponding to each point.

		Arguments:
		@@ -231,8 +265,8 @@ def LonLatPolygonCentroids( xx, yy, maxlevel=5, _level=0, indent='' ) :

		Return values:

		* ``xxp``, ``yyp`` : arrays of shape :math:`(P)` with the coordinates of the points
		* ``wwp`` : arrays of shape :math:`(P)` with area in m2 assigned to each point
		* ``xxp``, ``yyp`` : arrays of shape :math:`(P,N)` with the coordinates of the points
		* ``wwp`` : arrays of shape :math:`(P,N)` with area in m2 assigned to each point

		The polygon is recursively devided in triangles, the points are their centroids.
		First the centroid of the polygon is determined, and the first set of triangles is formed
		@@ -244,7 +278,8 @@ def LonLatPolygonCentroids( xx, yy, maxlevel=5, _level=0, indent='' ) :
		\|/ \\|
		o---o

		Each triangle is devided into 2 new triangles::
		Each triangle is devided into 2 new triangles with the middle of the longest side as their top
		and the two shortest sides as their base::

		2
		o
		@@ -301,9 +336,9 @@ def LonLatPolygonCentroids( xx, yy, maxlevel=5, _level=0, indent='' ) :
		#endif # triangle or polygon

		# return center:
		xxc = xc.copy()
		yyc = yc.copy()
		wwc = A.copy()
		xxc = xc.reshape((1,np))
		yyc = yc.reshape((1,np))
		wwc = A.reshape((1,np))

		else :

		@@ -325,9 +360,9 @@ def LonLatPolygonCentroids( xx, yy, maxlevel=5, _level=0, indent='' ) :
		if i == 0 :
		xxc,yyc,wwc = xxt,yyt,wwt
		else :
		xxc = numpy.hstack((xxc,xxt))
		yyc = numpy.hstack((yyc,yyt))
		wwc = numpy.hstack((wwc,wwt))
		xxc = numpy.vstack((xxc,xxt))
		yyc = numpy.vstack((yyc,yyt))
		wwc = numpy.vstack((wwc,wwt))
		#endif
		#endfor # corners

		@@ -372,19 +407,19 @@ def LonLatPolygonCentroids( xx, yy, maxlevel=5, _level=0, indent='' ) :
		xxt,yyt,wwt = LonLatPolygonCentroids( numpy.array([xxbase0,xxmid,xxtop]).T,
		numpy.array([yybase0,yymid,yytop]).T,
		_level=_level+1, maxlevel=maxlevel)
		# initialize result as copy:
		xxc = xxt.copy()
		yyc = yyt.copy()
		wwc = wwt.copy()
		# extend:
		xxc = xxt
		yyc = yyt
		wwc = wwt

		# ~ triangle 2
		xxt,yyt,wwt = LonLatPolygonCentroids( numpy.array([xxmid,xxbase1,xxtop]).T,
		numpy.array([yymid,yybase1,yytop]).T,
		_level=_level+1, maxlevel=maxlevel)
		# extend:
		xxc = numpy.hstack((xxc,xxt))
		yyc = numpy.hstack((yyc,yyt))
		wwc = numpy.hstack((wwc,wwt))
		xxc = numpy.vstack((xxc,xxt))
		yyc = numpy.vstack((yyc,yyt))
		wwc = numpy.vstack((wwc,wwt))

		#endif # multi or triangle