great circle arc spherical area about nco HOT 26 CLOSED

Done. This is in 4.5.2-beta1. Missing areas should be diagnosed automatically, no switch needed.
Need feedback on these test results: Areas of random gridcell (ncol=0) from native grid output and NCO diagnosis are 4.39677373284e-05 and 4.39677373279e-05 sr, respectively. Relative agreement to tolerance of 1.0e-10. NCO-diagnosed ne30 global area exactly equals original (native grid) ne30 global area = 12.5663706144 sr. These global areas equal 4_pi = 12.5663706128 sr to fractional tolerance of 1.0e-09. This seems low to me. Expected 4_pi agreement to about 1.0e-13 not 1.0e-09.
Implemented checks for ill-conditioned triangles (e.g., a=b=~0.5c) and found none.
Why do native grid cell vertices produce such poor areas? My understanding is that an offline program determines FE native grid vertices by nudging vertices around subject to global constraint on area. @mt5555 is this an artifact of the convergence criteria used to locate vertices on the FE grid? If so, can you re-run program and tighten global constraint to get closer to 4*pi?

from nco.

Which SCRIP file are you testing with? For our cubed-sphere grids, I think it should be better. For the CONUS grid, the version in the ACME inputdata only has an approximate dual grid.

from nco.

all tests done with ne30np4_pentagons.091226.nc
have not yet tried ne120, but would use ne120np4_pentagons.100310.nc.
have not yet tried CONUS.

from nco.

Nevermind. My "true value" for pi in the above comes is incorrect! True value is 12.56637061435917295376, and this agrees with native grid and NCO diagnosis. Will recompute tolerances based on true true value...

from nco.

Tolerances of global area from original ne30 pentagons SCRIP and NCO diagnosis from those vertices are 1.00000000000001082621 and 0.99999999999993761494, respectively. So both agree at 1.0e-13, as expected. Sorry for the false alarm! Are there are any more low-hanging regridding fruit that NCO can pick? If so, please say. Otherwise close this issue if it looks done to you.

from nco.

Recap: Yesterday I mentioned (correctly) that gridpoint area agreement between NCO and native grid area was within fractional relative tolerance of 1.0e-10, and that global integral agreement was similar. Turned out that I had used wrong value for 4pi and global integral agreement was within 1.0e-13 so we closed the ticket. I re-opened the ticket because I do not understand why the gridpoint agreement is so poor, i.e., only 1.0e-10. @mt5555 are you confident about the areas in the pentagons grid file?
´´´
zender@roulee:~$ ncks -H -s %20.15e -v grid_area -d grid_size,0 ${DATA}/grids/ne30np4_pentagons.091226.nc
4.396773732840135e-05

´´´
NCO gets 4.396773732788860e-05 for this same area (the first ne30 gridcell). Why the difference? Is this the great circle arc issue? NCO uses great circle arcs and I thought the native grid made the same assumption. Does it?

from nco.

The SCRIP file for the HOMME spectral element grid was generated internally by HOMME and uses Girad's formula to compute areas. So for an area of size 4.e-5, the area computed by HOMME could loose up to 5 digits of accuracy. Hence my first guess is the error is in the SCRIP file produced by HOMME.

There's a newer grid file, produced by our Matlab code. Can you try:

https://acme-svn2.ornl.gov/acme-repo/acme/mapping/grids/conusx4v1np4_chevrons_scrip_c150815.nc

from nco.

For ${DATA}/grids/conusx4v1np4_chevrons_scrip_c150815.nc, NCO-diagnosed and raw grid area of first gridcell:
zender@roulee:$ ncks -H -s %20.15e -v area -d ncol,0 ${DATA}/ne30/rgr/dogfood.nc
3.653857995289623e-05
zender@roulee:$ ncks -H -s %20.15e -v grid_area -d grid_size,0 ${DATA}/grids/conusx4v1np4_chevrons_scrip_c150815.nc
3.653857995295246e-05

Ratio is 3.653857995289623/3.653857995295246=0.99999999999846107867, i.e., agreement to fractional relative tolerance of 1.0e-11, almost 1.0e-12. Better than the ne30 grid previously tested, yet underwhelming. Would expect 1.0e-14 agreement or better if NCO and the Matlab script both
implement L'Huillier's formula corrrectly with same inputs yet different order of operations. Have any hypotheses consistent with these tests (aside from someone goofed)?

from nco.

I'll get that same number output from the matlab code.

from nco.

3.653857995294302e-05 area computed via matlab routine
3.653857995295246e-05 (GL weight)
relative error between the above two numbers: 3e-13.
Relative error between NCKS and Matlab area calculation: 1.3e-12.

For the record: the matlab routine searches for a dual grid with spherical area that matches the GL weight. It uses a Newton iteration, and if we run more iterations it wont converge to better than 3e-13.

The matlab area is computed via a convoluted process (that could make it less accurate than ncks)
0. assume the SCRIP control volume is a "n-gon"

1. convert lat/lon coordinates to cartesian
2. compute a mean of all the coordinates (centroid of the n-gon)
3. using the centroid, and each of the edges, make n triangles
4. sum the area of the n triangles (via L'Huillier's formula)

So the matlab code is converting to Cartesian, and then in L'Huillier's formula, converting back to lat/lon to compute great circle distances). And our matlab algorithm is also forming n triangles, and ncks might make triangles out of just vertices and you'd only have n-1 triangles?

This brings up another possibility: some of the cells in our dual grid are (slightly) not convex, and ncks may be assuming they are convex when breaking up into triangles in order to compute the spherical area? That's why we switched to using a centroid, since it's a little more robust.

from nco.

Making progress...there seem to be enough differences to explain the discrepancies.

To summarize, my understanding is that:
Matlab starts with the GL weight and shifts vertices around to match.
Computes area based on centroid and N triangles.
Also the "areas" listed in the SCRIP file are actually GL weights, and the vertices in the SCRIP file are the polygons that best fit those weights (and are known to be hard/impossible to improve).
Yes?

NCO assumes the vertices in the SCRIP file are gospel and computes (via L'Huillier) the area using N-2 (not N-1) triangles.
Is L'Huillier applied to N centroid-based triangles more accurate than applied to N-2 vertice-based triangles?
Because vertice-based triangles are "sharper" than centroid-based?
Is this known "for sure" or maybe a coincidence for this point?
Surpised that using fewer triangles appears to be less accurate!
Would naively expect less round-off noise from fewer triangles.

What, if anything should be done differently?
NCO could compute a centroid and use N triangles, like Matlab.
Or use N-2 triangles with a side-angle-side formula.
Or keep N-2 triangles and L'Huillier.
However, if the centroid algorithm is known to be more accurate, then NCO should just use that.

from nco.

Yes to the first question (areas in the file are actually the GL weights).

I would think the ncks formula (using N-2 triangles) is more accurate. I could switch to it in the matlab code, if you have a way to handle non-convex polygons? i.e. consider a chevron. Some of the N-2 triangles may not be inside the chevron?

from nco.

from nco.

Good point. NCO algorithm with N-2 triangles only works for convex polygons. Not an issue for the point at ncol=0 (only four vertices), but definitely an issue elsewhere in that grid. Ran out of time today...will consider options during REM.

from nco.

The NCO algorithm (i.e., N-2 triangles) always works for N=3 (obviously), and for N>=4 convex polygons. The NCO algorithm can be made to work for N=4 concave polygons by starting at the vertex that contains the concave angle. The NCO algorithm would break-down into a morass of if/thens for concave polygons N>=5.

The centroid N-triangle algorithm always works for N<=4, for concave and convex, and always for N>=5 for convex. It breaks down for (some) concave polygons N>=5. It does work for the N=6 chevron you drew. What other concave shapes are used/important?

Correct me if I'm wrong but we expect fewer triangles (i.e., N-2) to be (slightly) more accurate than more triangles because of less round-off error. (We still do not understand why NCO and matlab produce such different results for ncol=0 on the conus grid).

The general case for N>=5 requires lots of trig, for both algorithms. I'm not interested in NCO solving the general case of polygonal area. I am interested in NCO solving polygonal area for all shapes that arise in climate modeling grids used by funding agencies :)

One way NCO could proceed is to require a "check_concave" switch. When not activated, NCO would use its current fast (N-2 triangle, no need to solve for interior angles), accurate (less round-off error) method. When activated, NCO would check for concave angles for all N>=4 polygons. This requires solving all interior angles (whereas for convex polygons NCO only solves for side/arc lengths). If polygon is convex, use L'Huillier N-2 triangle solution. If polygon is concave, use N-triangle centroid solution. Is Girard's formula valid for concave polygons? If so that is also an option, though it sacrifices precision. Yet if Girard works in the general case, then it might be "good enough" if restricted to the (presumably very few) concave polygons on an otherwise predominantly convex grid.

I'm most interested in getting this right for the Cubed Sphere (CS) grid at the moment, with RRM grids (and concave shapes like chevrons etc.) to follow. My understanding is that all shapes on the default CS grid (i.e., default dual grid for CAM SE) are convex, mainly quadrilaterals, with maybe pentagons/hexagons at the edges. Is this correct?

from nco.

lets track down the loss of precision for this first cell and see if it is related to using N or N-2 triangles. I will compute both ways in our matlab code and see if one agrees with NCO. If they two different calculations in Matlab are different, I would assume the N-2 version is more accurate.

In all cases, I vote for L'Huillier's formula. Both Girard and L'Huillier are for triangles, so both require decomposing into triangles and wont fix the non-convext issue. Girard's formula will compute (area+pi), and then subtract pi. So if area is 1e-5, subtracting pi will loose 5 digits of accuracy.

For the record: a spectral element grid has all quad elements - but it's dual grid is not well defined and our simplest algorithm to compute the dual grid will make chevrons that are only slightly convex. So I think its a fair assumption to assume the center of mass is always inside the cell. What about (N-2) algorithm for N<=4, and the N algorithm for N>4?

from nco.

Thank you for doing this. I should implement the centroid rule as an option in NCO and verify it agrees with the Matlab answer. By "Girard's formula" I'm referring to the extended version that applies to N-gons, not just to triangles. I don't know what it's real name is, but it looks just like an extended version of Girard's formula for N-gons. See bottom of Wikipedia page on Spherical Triangles. I suspect that it only applies to convex polygons.

Re: accuracy, the supplied vertices on the dual grid are only accurate to ~13 decimal places. Losing three more digits at possibly only a few points may be a worthwhile tradeoff to solve the general case. Local results lose precision, yet global integrals would not change much. Anyway, I suggest considering it if a) Generalized Girard applies to all N-gons including concave, and b) there are grids containing weird shapes for which the centroid algorithm does not work (I assume you meant "slightly concave" above).

N-2 algorithm for N<=4 and N-algorithm for N>4 seems like a good compromise. Does the dual-grid ever have concave quadrilaterals? If so, need to revise N-2 algorithm to start with the right apex point. I will implement N-algorithm for N>4 and independently see if it reproduces matlab results.

from nco.

The N and N-2 formulas appear equally accurate:
GLL = 3.653857995295246e-05
N-2 = 3.653857995294302e-05
N triangles = 3.653857995294301e-05
nco: 3.653857995289623e-05

Both the matlab code and the NCO code are reading in the same coordinates from the netcdf file, so I think the difference has to be in our implementation of L'H's formula, or in the trig functions used by matlab vs. the NCO C compiler?

Based on this result (and for simplicity), I would argue for using the centroid algorithm for all polygons, unless it impacts the speed?

from nco.

just for the record:
lat1 =
0.341215693632680
0.341147888120055
0.346916154821447
0.347227105822592
0.347227105822592
0.347227105822592

lon1 =
3.815996778621659
3.822830543733946
3.819605294803501
3.813060805743338
3.813060805743338
3.813060805743338

matlab code:

a = grtcircdist_latlon(lat1(1),lon1(1),lat1(2),lon1(2));
b = grtcircdist_latlon(lat1(2),lon1(2),lat1(3),lon1(3));
c = grtcircdist_latlon(lat1(3),lon1(3),lat1(1),lon1(1));
s = (a+b+c)/2;
A1 = sum(4_atan(sqrt(tan(s/2)._tan((s-a)/2)._tan((s-b)/2)._tan((s-c)/2))));

a = grtcircdist_latlon(lat1(1),lon1(1),lat1(3),lon1(3));
b = grtcircdist_latlon(lat1(4),lon1(4),lat1(3),lon1(3));
c = grtcircdist_latlon(lat1(4),lon1(4),lat1(1),lon1(1));

s = (a+b+c)/2;
A2 = 4_atan(sqrt(tan(s/2)._tan((s-a)/2)._tan((s-b)/2)._tan((s-c)/ ...
2)));

A1+A2

from nco.

Regrading the wikapedia formula for the N-gon: I wouldn't trust it, since it is of the form:

A1 = (N-2)*pi
A2 = A1 + area

and one computes A1 and A2 and then subtracts them, loosing precision when area<<1

from nco.

function sigma = grtcircdist_latlon(lat1,lon1,lat2,lon2)

delphi = abs(lat1 - lat2);
dellam = abs(lon1 - lon2);

sigma = 2_asin(sqrt((sin(delphi/2)).^2 + cos(lat1)._cos(lat2).*(sin(dellam/2)).^2));

from nco.

    cos_a=lat_bnd_cos[idx_a]*lon_bnd_cos[idx_a]*lat_bnd_cos[idx_b]*lon_bnd_cos[idx_b]+
      lat_bnd_cos[idx_a]*lon_bnd_sin[idx_a]*lat_bnd_cos[idx_b]*lon_bnd_sin[idx_b]+
      lat_bnd_sin[idx_a]*lat_bnd_sin[idx_b];
        arc_a=acos(cos_a);

from nco.

Take a look at:

https://en.wikipedia.org/wiki/Great-circle_distance

They seem to be saying the arcsin formula used in our Matlab code is more accurate than the arccos formula used by NCO (cut & paste below).

On computer systems with low floating-point precision, this formula (the arccos formula) can have large rounding errors if the distance is small (if the two points are a kilometer apart on the surface of the Earth, the cosine of the central angle comes out 0.99999999). For modern 64-bit floating-point numbers, the spherical law of cosines formula, given above, does not have serious rounding errors for distances larger than a few meters on the surface of the Earth.[2] The haversine formula (the arcsin formula) is numerically better-conditioned for small distances:[3]

from nco.

Yes, that explains most of the discrepancy. Implemented haversine formula for great circle arc lengths and tested against ncol=0 on conus chevrons file (same as before):
3.653857995295246e-05 raw GLL weight
3.653857995294302e-05 matlab N-2 triangles
3.653857995294301e-05 matlab N triangles
3.653857995294258e-05 new NCO (haversine)
3.653857995289623e-05 old NCO (acos)
Still a small discrepancy in favor of matlab code. Possibly order of operations? Will check if duplicate same order as matlab...

from nco.

Nope, order of operations does not explain remaining discrepancy.

from nco.

relative error between NCO and Matlab implementations: 1.2e-14

That seems pretty reasonable...

We could repeat the calculation in quad precision and determine which formula was more accurate, but I'm not sure how the less accurate one could be improved?

from nco.