CubedSphere grid cell area calculated by exchange grid tool has discontinuous values around the two poles

[nikizadehgfdl]

Describe the bug
The cubed sphere grid cell area calculated by the exchange grid tool has bumps around the two poles (see figures below). These bumps are not there in the cell area calculated by make_hgrid itself, the tool that creates the CS grid. Note that make_hgrid uses spherical_excess_area() (Gauss-Bonnet formula area_of_polygon = sum of interior spherical angles - 2pi) , whereas the exchange grid tool uses poly_area() ( area_of_polygon = line_intergal sin(lat) d(lon) ). The question is which one of these formulas are more accurate near the singular poles?

Non-GCA Stretched grids have the same symptoms which results in tiling errors.

To Reproduce
Get a version of frenctools with more verbosity for the exchange grid creation tool:

bash
git clone -b nnz_NONGCA_strechedGridIssues https://gitlab.gfdl.noaa.gov/fre/fre-nctools.git fre-nctools_gitlab
cd fre-nctools_gitlab
git checkout adde039  #adds calculated atm_area to land_mask files
export CONFIG_SITE=site-configs/gfdl-ws/config.site 
source site-configs/gfdl-ws/env.sh 
module load gcc/6.2.0
export PATH=.:$PATH
autoreconf -i
./configure 
make

Now create any coupled mosaic grid with a uniform CS grid, e.g., Test03-grid_coupled_model.sh
Plot the calculated area_atm in land_mask_tile6.nc

Expected behavior
computed area in C48_grid.tile6.nc does not have a bump at the pole.

System Environment
FRE-NCTOOLS build on PAN

Additional context
We already know that the poly_area() formula is not correct when side of polygon crosses a pole.

Figure below shows area of ATM grid cells along the longitude passing through pole (at the peak) calculated by the exchange grid tool (poly_area) for c256_om4p25 grid

Same ATM grid cell area as above calculated by the make_hgrid.c (spherical_excess_area())

[ngs333]

@nikizadehgfdl
I am not sure that this observation can help some users, but while looking at the poly_area() discontinuities, I noticed that make_coupler_mosaic will not create the singularities if make_hgrid uses the gca for all the grids
and make_coupler_mosaic is given the argument "--interp_order 1". Indeed, the discontinuity is observed in file "land_mask_tile6.nc" created by running Test03. However, if we change Test03 so that
" --great_circle_algorithm" is provided as an argument to (all calls of) make_hgrid, and "--interp_order 1" is provided as argument to make_coupler_mosaic, the ncview plots of the central point of field "area_lnd" changes from the plot on the left to the one on the right in the attached figure:

[nikizadehgfdl]

@ngs333 thanks for the observation. As you noted, the problem with GC algorithm is that currently it only works with 1st order interpolation. Whereas almost everyone wants to use 2nd order conservative method. Also, it produces non-convex exchange cells for stretched grids which renders it totally useless with our flux-exchange code (models crash).

I suspect that the bump in the cell area is due to a bug and is not a feature of these grids because the exact Gauss-Bonnet formula used in make_hgrid.c does not show it. Perhaps when you use GC algorithm the exact formula is used (rather than a call to the function poly_area())? I'll probably get to this next week.

[ngs333]

@nikizadehgfdl
Thanks! Also, do you have at this point any further intuition as to the awnser to your question : " ... which one of these formulas ( spherical_excess() or poly_are() )are more accurate near the singular poles? " ?

I'll try to see if I can find a bug in poly_area(). But I am also thinking of changing make_coupler_mosaic to be able to use spherical_excess() in lieu of poly_area(). Do you think that could lead to something usefull?

[nikizadehgfdl]

I think the spherical_excess() is the accurate geometrical formula everywhere on sphere (discovered by Lambert in 1760's for sphere and extended by Gauss and Bonnet to arbitrary surfaces). I don't know why poly_area() was created/used instead of spherical_excess(). I think your suggestion of using the spherical_excess() in make_coupler_mosaic is in the right direction.

[ngs333]

Update from week of Jan 3, 2022:
Intel Inspector showed no issues (like invalid_memory access or uninitalized memory access) in make_coupler_mosaic that can cause wrong awnsers. Also, I've spend some time stepping through
poly_area() as called my make_coupler_mosaic and found no issues, bu this will continue.

[ngs333]

@nikizadehgfdl
Niki, Thanks for your email response today. Also, I have some questions (which may be better covered in a meeting) concerning function poly_area.

In email you clarified
"I think the excess angle formula (Lambert: area = (angle sum-pi)*R**2 ) only works if the sides of the triangle are geodesics (great circle sections).
If the sides are arbitrary curves (as might be the case for exchange grid cells) we need to integrate."

Might we need a third function as poly area has too many assumptions for its use with the exchange grids at issue? Consider the signature of the function:
double poly_area(const double x[], const double y[], int n) .
x[] and y[] are in radias, and it assumed a fixed R! Would this be too restrictive for use for the exchange grid cells in question? Could t this ( https://math.stackexchange.com/questions/3207981/caculate-area-of-polygon-in-3d ) be a candidate 3rd formula ?
Thanks

[ngs333]

@nikizadehgfdl
Looks like the singularity of function poly_area() is not do to the grids being stretched but merely do to the cells having a pole point within them. I should have noticed this before when I ran some of the tests I mention above. However, there is some related literature from recent years that can help, plus also I have a suggestion I am trying to do in “pencil and paper” and should lend to as easy algorithm.

First there is the publication “Stokes’s Theorem, Data, and the Polar Ice Caps”, by Baryshnikov and Ghrist that hose some new formula (see formula 6) that should be super accurate, plus they outline how to deal with the pole singularities (see last section). Finally, they mention another, slightly older publication that may also be of help:

Chamberlain, R. G., Duquette, W. H. (2007). Some algorithms for polygons on a sphere. JPL Publication 07-3. Pasadena, CA: JPL NASA.

Reading the Baryshnikov paper gave me an idea, which is when a polygon encloses a pole, calculate its are by dividing the polygon in two and summing the areas of the two but in the manner the image describes . In the attached PDF, the two polygons (one is a triangle) are formed by inserting two additional lines : one between the pole and any point (say P), and a second between the pole and the the point preceding (of following) P in circular order. For those two new sub-polygons, function poly_area can be used. Or perhaps
another similar function if not poly area.
split2.pdf

[nikizadehgfdl]

Yes, this bug has nothing to do with stretched grid, it was seen in regular CS c256 grid.It's probably a bug in poly_area() formula when an edge of the grid cell passes through a pole. I am not sure if the same discontinuity happens when grid cell totally encloses a pole (as your figures on the left show).

[ngs333]

@nikizadehgfdl, @bensonr , @ceblanton , @jwdGFDL
Hi, based upon some comments from Jeff in my 2022 work plan, I am adding this rough outline of possible steps
on this issue :
1A) Identify all conditions where area is calculated incorrectly. (e.g. area is enclosing pole, are near poles but not enclosing, both, other areas)
1B) Produce a file with sample polygon for each condition above.
1C) For test purposes, are there any polygons/triangles to add to the list because their areas are known exactly or analytically?
2A) Inspect and compare the area calculations made by make_hgrid (MH) and make_coupler_mosaic (MCM) for each of those polygons in 1B).
2B) Is MH and MCM really performing the poly_area calculations on the exact same polygons and with the exact same function (poly_area) in tests (e.g. Test03 ).
3A) Run Intel_inspector on MCM
3B) Run MCM in debugger and look at the calculation of erroneous areas.
3C) Identify how poly_area can be improved.
4A) Identify alternatives to formula for poly area (See especially formulas in paper by Baryshnikov, in paper by Chamberlain, projecting onto plane via area-preserving projection, splitting the rectangle into triangles (see above diagrams). Code and compare results to poly_area.
4B) Identify other packages (e.g. python packages) and see how well they calculate the area for the sample polygons.
5A) Depending on how well poly_area and MCM is fixed before step 4, and how accurately other formulas calculate the areas, consider creating another NCTools github issue/enhancement for a 3rd method of polygon area calculations NCTools.
6A) Make unit tests to make sure the stretched grids workflow are producing accurate results for polygon areas.
6B) Make unit tests that compare the polygon areas of MH vs those of MCM.
6C) Make utility that compare the polygon areas of MH vs those of MCM.

MZ Current (1/19/2022) status: I am starting again from 1A. I had done 1A, 2A, 3A, 3B, 4C but I think now that I must have made some mistakes or missed something.

[ngs333]

@nikizadehgfdl, @bensonr , @ceblanton , @jwdGFDL

Is make_coupler_mosaic making "five sided squares" at the poles?
It seems most squares are specified by 4 points, like
latitudes : [ -35.2644, -35.9889, -36.7609, -35.9889 ]
longitudes : [215, 213.427, 215, 216.573]
and they appear as 4-sided rectangles as in the image below:

But at a pole they are being specified by 5 points :
lattitudes: [ -90.,-90.,-87.9225, -87.0632, -87.9225]
longitudes :[-10., 80., 80., 35., -10.]
and appear as regular squares as in the image below:

An extra line in the path of the integral can be the extra value that the area is getting.

[ngs333]

Some data and notes :

1. For unit test Test03, File poly_5p_sides_all_v2.txt (attached) has the 20 polygons procesed by make_coupler_mosaic that are five sided (though they appear 4 sided as in the figure above).
2. The five sided polygons are made by function fix_lon() that prepares polygons for use with poly_area() so that
   poly_area() can get the right awnser. (This also occurs in other apps, like make_hgrid). It begins by checking the latitude of each polygon point, and if the latitude is too close to +-90, it modiffues the polygon. In all cases it turned a four-sided polygon into a five-sided polygon. It is suspicious that the resultant polygon still has some latitudes equal to +-90. It might be worth investigating the making of a different polygon, or slightly altering a polygon point differently. For example, such as setting the lat of the point with lat == 90, to lat = (90 - epsilon). Everything else of the polygon can be kept the same. The calcualted area may be smaller by about a factor relating to epsilon, but the final area can be adjusted and may be much better than what is currently there.
   poly_5p_sides_all_v2.txt

[ngs333]

@nikizadehgfdl
I have looked at poly_area() formula for a bug but I have not found it. I have also tried changing fix_lon() (such as by
making a different but either similar or equivalent polygon) but this has not helped the calculation in poly_are. One thing that I have tried that may be worth mentioning though it does't seem to work quite "well enough" is to circumvent both fix_lon() and
poly_area() when (only when) a polygon has a point at a pole (i.e. with lat being equal to +- 90 within epsion). In these cases the spherical excess area is used. The result is that instead of the discontinuity showing plot above ( at #44 (comment) ) we get this plot below:

The plot does not show the discontinuity, but there seems to be a dip at the pole.

[ngs333]

@nikizadehgfdl
I added exploratory code (https://github.com/ngs333/FRE-NCtools/tree/polyarea1) that rotates polygons that have a
a point at a pole (i.e. with lat being equal to +- 90 within epsilon). More specifically, when a polygon area is needed, if the polygon is close to the pole, a copy of the polygon is made and it is rotated about a pre-determined axis
so that it is away from the pole. When poly_area() is applied to the rotated polygons (and fix_lon() is not applied),
the calculated area , for a targeted observed set, is roughly 0.1% less than what the spherical excess area function
calculated. ( For the points of discontinuity on the plots way above, that original fix_lon() + poly_area() calculates
areas that are about 20% greater than the spherical excess area.)

[ngs333]

I looks like rotating (away from the poles) the two polygons that are input arguments to the polygon clipping algorithm
should be tried as a potenntiall fix to the remaining bugs. The figure on the top is for a clipping
of the common case where the two polygons are away from the poles. The input polygons are yellow and red, the resultant clip is in purple with dark dashed lines for the edge. From inspection it is obviously a good clip. The other figure shows a clip
wiith near-pole input polygons. The resultant clip is not the intersection of the input polygons, and it also seems to not have the right shape and area.

(PS The graphics code is in Python using numpy and matplot3d. It is in my gfdl_misc repository, directory src/python, file polydis3d.py (or https://github.com/ngs333/gfdl_misc/blob/master/src/python/polydis3D.py )

[nikizadehgfdl]

@ngs333, this is odd. What are the two grids you are using and what are the coordinates of the corners of the two cells? Does this happen with the bronx-19 fre-nctools exec?

[ngs333]

@nikizadehgfdl
Yes it is odd. I don't know right now if this happens with bronx-19. In this case, the code base has been
modified to calculate areas of a polygon near a pole by rotating a copy of the polygon away from the pole. Function fix_lon is correspondingly modified by not exercising that part that adds and/or removes vertices - i.e. it really only fixes the lons. So in principle the mods to the fix_lon could have changed the inputpolygons to the clipping functions. In any case, I am suspicious of the clipping functions and also fix_lon in original form. The input polygons are from the ATMxLND and OCN. Thats Lines 1617 ff in the current baseline (https://github.com/NOAA-GFDL/FRE-NCtools/blob/master/tools/make_coupler_mosaic/make_coupler_mosaic.c ) with the print functions modified to print the xyz (not the lat/lon) coordiantes with earh radius = 1. I can work on this again on Wed and do other tests if you like. Thanks.

[ngs333]

@bensonr , @nikizadehgfdl
I tried the idea of rotating the two input argument polygons of clip_2dx2d away from the poles whenever any point of either
of the two polygons was within 0.5 deg from a pole. It got rid of a lot of problems but there is still something wrong somewhere.

This is the current tail end of the make_coupler_mosaic report:

The maximum area mismatch is at tile=6, i=24, j=25, ratio=0.0995906, area_atm_grid=5.36158e+10, area_atm_xgrid=4.82762e+10
NOTE: axo_area_sum is 378520935785634.187500 and ocean fraction is 74.220077%
NOTE: axl_area_sum is 131477002189518.593750 and land  fraction is 25.779904%
NOTE: tiling error is 0.000019%

The tiling error above is more than two orders of magnitude better than without the rotation.

The most interesting thing is that tile 3 is now almost completely devoid of errors, while tile 6 is a little worse!
Before they were similarly bad.

[ngs333]

This is just to note the effect of the changes (latest code is now here : https://github.com/ngs333/FRE-NCtools/tree/polyarea2 ) with the C256 streched grid Niki suggested. The changes include:
a) polygons with any vertex closer than 0.5% from a pole are rotated away from the pole.
b) Function fix-lon is modified to only change the longitudes (so the addition and removal of vertices to make function
poly_area work near poles is not done)
Although the polygon area discontinuity at the poles are not present (i.e. polygons touching the poles dont have the huge areas anymore) , make_coupler_mosaic shows about a 1% discrepancy bewteen areas at the pole via comparison of the AMT grid cell area at the pole vs the sum of exchange grid cell areas clipped with same ATM grid cell area. The tail of of the output from make_coupler_mosaic is the following

at tile =6, i=127, j=130, ratio=0.00875916, area1=1.87636e+09, area2=1.8928e+09
at tile =6, i=128, j=130, ratio=0.00823054, area1=1.87636e+09, area2=1.89181e+09
at tile =6, i=129, j=130, ratio=0.00523273, area1=1.87619e+0 9, area2=1.86637e+09
The maximum area mismatch is at tile=3, i=128, j=127, ratio=0.10891, area_atm_grid=1.87688e+09, area_atm_xgrid=1.67247e+09

NOTE: axo_area_sum is 360506610213201.562500 and ocean fraction is 70.678971%
NOTE: axl_area_sum is 149553928650354.843750 and land  fraction is 29.320732%
NOTE: tiling error is 0.000297%