1 /**
2 * @file area.c
3 *
4 * @copyright Copyright (C) 2013 Moritz Hanke <hanke@dkrz.de>
5 * Rene Redler <rene.redler@mpimet.mpg.de>
6 *
7 * @version 1.0
8 * @author Moritz Hanke <hanke@dkrz.de>
9 * Rene Redler <rene.redler@mpimet.mpg.de>
10 */
11 /*
12 * Keywords:
13 * Maintainer: Moritz Hanke <hanke@dkrz.de>
14 * Rene Redler <rene.redler@mpimet.mpg.de>
15 * URL: https://dkrz-sw.gitlab-pages.dkrz.de/yac/
16 *
17 * This file is part of YAC.
18 *
19 * Redistribution and use in source and binary forms, with or without
20 * modification, are permitted provided that the following conditions are
21 * met:
22 *
23 * Redistributions of source code must retain the above copyright notice,
24 * this list of conditions and the following disclaimer.
25 *
26 * Redistributions in binary form must reproduce the above copyright
27 * notice, this list of conditions and the following disclaimer in the
28 * documentation and/or other materials provided with the distribution.
29 *
30 * Neither the name of the DKRZ GmbH nor the names of its contributors
31 * may be used to endorse or promote products derived from this software
32 * without specific prior written permission.
33 *
34 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
35 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
36 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
37 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
38 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
39 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
40 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
41 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
42 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
43 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
44 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
45 */
46
47 #include <stdlib.h>
48 #include <stdio.h>
49 #include <math.h>
50 #include <float.h>
51
52 #include "area.h"
53 #include "grid.h"
54 #include "clipping.h"
55 #include "geometry.h"
56 #include "utils.h"
57 #include "ensure_array_size.h"
58
59 static inline double scalar_product(double a[], double b[]);
60
61 static inline double inner_angle ( double plat, double plon,
62 double qlon, double qlat );
63
64 /* ----------------------------------- */
65
yac_triangle_area(struct grid_cell cell)66 double yac_triangle_area ( struct grid_cell cell ) {
67
68 /* taken from the ICON code, mo_base_geometry.f90
69 provided by Luis Kornblueh, MPI-M. */
70
71 const double tol = 1e-18;
72
73 double s01, s12, s20;
74 double ca1, ca2, ca3;
75 double a1, a2, a3;
76
77 double * triangle[3];
78 double u01[3], u12[3], u20[3];
79
80 if ( cell.num_corners != 3 ) {
81 printf ("Only for triangles!\n");
82 return -1;
83 }
84
85 triangle[0] = cell.coordinates_xyz[0];
86 triangle[1] = cell.coordinates_xyz[1];
87 triangle[2] = cell.coordinates_xyz[2];
88
89 /* First, compute cross products Uij = Vi x Vj. */
90
91 crossproduct_ld(triangle[0], triangle[1], u01);
92 crossproduct_ld(triangle[1], triangle[2], u12);
93 crossproduct_ld(triangle[2], triangle[0], u20);
94
95 /* Normalize Uij to unit vectors. */
96
97 s01 = scalar_product(u01, u01);
98 s12 = scalar_product(u12, u12);
99 s20 = scalar_product(u20, u20);
100
101 /* Test for a degenerated triangle associated with collinear vertices. */
102
103 if ( fabs(s01) < tol &&
104 fabs(s01) < tol &&
105 fabs(s01) < tol )
106 return 0.0;
107
108 s01 = sqrt(s01);
109 s12 = sqrt(s12);
110 s20 = sqrt(s20);
111
112 for (int m = 0; m < 3; m++ ) {
113 u01[m] = u01[m]/s01;
114 u12[m] = u12[m]/s12;
115 u20[m] = u20[m]/s20;
116 }
117
118 /* Compute interior angles Ai as the dihedral angles between planes:
119
120 CA1 = cos(A1) = -<U01,U20>
121 CA2 = cos(A2) = -<U12,U01>
122 CA3 = cos(A3) = -<U20,U12>
123
124 */
125
126 ca1 = -u01[0]*u20[0]-u01[1]*u20[1]-u01[2]*u20[2];
127 ca2 = -u12[0]*u01[0]-u12[1]*u01[1]-u12[2]*u01[2];
128 ca3 = -u20[0]*u12[0]-u20[1]*u12[1]-u20[2]*u12[2];
129
130 if ( ca1 < -1.0 ) ca1 = -1.0;
131 if ( ca1 > 1.0 ) ca1 = 1.0;
132 if ( ca2 < -1.0 ) ca2 = -1.0;
133 if ( ca2 > 1.0 ) ca2 = 1.0;
134 if ( ca3 < -1.0 ) ca3 = -1.0;
135 if ( ca3 > 1.0 ) ca3 = 1.0;
136
137 a1 = acos(ca1);
138 a2 = acos(ca2);
139 a3 = acos(ca3);
140
141 /* Compute areas = a1 + a2 + a3 - pi.
142
143 here for a unit sphere: */
144
145 // return MAX(( a1+a2+a3-M_PI ) * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS, 0.0);
146 return MAX( (a1+a2+a3-M_PI) , 0.0 );
147 }
148
149 /* ----------------------------------- */
150
yac_cell_area(struct grid_cell cell)151 double yac_cell_area ( struct grid_cell cell ) {
152
153 /* generalised version based on the ICON code, mo_base_geometry.f90
154 provided by Luis Kornblueh, MPI-M. */
155
156 const double tol = 1e-18;
157
158 size_t const M = cell.num_corners; // number of vertices
159
160 double area, s[M], ca[M], a[M], * p[M], u[M][3];
161
162 for (size_t i = 0; i < M; ++i) p[i] = cell.coordinates_xyz[i];
163
164 /* First, compute cross products Uij = Vi x Vj. */
165
166 for (size_t m = 0; m < M; m++ ) crossproduct_ld (p[m], p[(m+1)%M], u[m]);
167
168 /* Normalize Uij to unit vectors. */
169
170 area = 0.0;
171
172 for (size_t m = 0; m < M; m++ ) {
173 s[m] = scalar_product(u[m], u[m]);
174 area += s[m];
175 }
176
177 /* Test for a degenerated cells associated with collinear vertices. */
178
179 if ( fabs(area) < tol ) {
180
181 for (size_t m = 0; m < M; m++ ) s[m] = sqrt(s[m]);
182
183 for (size_t m = 0; m < M; m++ )
184 for (size_t i = 0; i < 3; i++ )
185 u[m][i] = u[m][i]/s[m];
186
187 /* Compute interior angles Ai as the dihedral angles between planes
188 by using the definition of the scalar product
189
190 ab = |a| |b| cos (phi)
191
192 As a and b are already normalised this reduces to
193
194 ab = cos (phi)
195
196 There is no explanation so far for the - in the loop below.
197 But otherwise we don't get the correct results for triangles
198 and cells. Must have something to do with the theorem.
199
200 */
201
202 for (size_t m = 0; m < M; m++ ) {
203 ca[m] = - scalar_product(u[m], u[(m+1)%M]);
204 if ( ca[m] < -1.0 ) ca[m] = -1.0;
205 if ( ca[m] > 1.0 ) ca[m] = 1.0;
206 a[m] = acos(ca[m]);
207 }
208
209 /* Compute areas = a1 + a2 + a3 - (M-2) * pi.
210
211 here for a unit sphere: */
212
213 area = - (double) (M-2) * M_PI;
214
215 for (size_t m = 0; m < M; m++ )
216 area += a[m];
217 }
218
219 // return MAX(area * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS, 0.0);
220 return MAX(area,0.0);
221 }
222
223 /* ----------------------------------- */
224
225 static inline struct sin_cos_angle
tri_area_quarter_angle(struct sin_cos_angle angle)226 tri_area_quarter_angle(struct sin_cos_angle angle) {
227
228 // the angle passed to this routine should be in the range [0;3*PI/4],
229 // in case the angle is outside this, we have to assume that this is due to
230 // numerical inaccuracy or
231 // tri_area was called for a too big triangle...has to be really big...
232
233 if (compare_angles(angle, SIN_COS_7_M_PI_4) >= 0) return SIN_COS_ZERO;
234 else return quarter_angle(angle);
235 }
236
237 /** area of a spherical triangle based on L'Huilier's Theorem
238 *
239 * source code is taken from code by Robert Oehmke of Earth System Modeling
240 * Framework (www.earthsystemmodeling.org)
241 *
242 * it has been extended by a more accurate computation of vector angles
243 *
244 * the license statement for this routine is as follows:
245 * Earth System Modeling Framework
246 * Copyright 2002-2013, University Corporation for Atmospheric Research,
247 * Massachusetts Institute of Technology, Geophysical Fluid Dynamics
248 * Laboratory, University of Michigan, National Centers for Environmental
249 * Prediction, Los Alamos National Laboratory, Argonne National Laboratory,
250 * NASA Goddard Space Flight Center.
251 * Licensed under the University of Illinois-NCSA License.
252 *
253 * \remark all edges are on great circle
254 */
tri_area(double u[3],double v[3],double w[3])255 static double tri_area(double u[3], double v[3], double w[3]) {
256
257 struct sin_cos_angle angle_a = get_vector_angle_2(u,v);
258 struct sin_cos_angle angle_b = get_vector_angle_2(u,w);
259 struct sin_cos_angle angle_c = get_vector_angle_2(w,v);
260
261 if (compare_angles(angle_a, SIN_COS_LOW_TOL) < 0) return 0.0;
262 if (compare_angles(angle_b, SIN_COS_LOW_TOL) < 0) return 0.0;
263 if (compare_angles(angle_c, SIN_COS_LOW_TOL) < 0) return 0.0;
264
265 double sin_sin = angle_a.sin * angle_b.sin;
266 double sin_cos = angle_a.sin * angle_b.cos;
267 double cos_sin = angle_a.cos * angle_b.sin;
268 double cos_cos = angle_a.cos * angle_b.cos;
269
270 double sin_sin_sin = sin_sin * angle_c.sin;
271 double sin_sin_cos = sin_sin * angle_c.cos;
272 double sin_cos_sin = sin_cos * angle_c.sin;
273 double sin_cos_cos = sin_cos * angle_c.cos;
274 double cos_sin_sin = cos_sin * angle_c.sin;
275 double cos_sin_cos = cos_sin * angle_c.cos;
276 double cos_cos_sin = cos_cos * angle_c.sin;
277 double cos_cos_cos = cos_cos * angle_c.cos;
278
279 double t_sin_a = sin_sin_sin - sin_cos_cos;
280 double t_sin_b = cos_sin_cos + cos_cos_sin;
281 double t_sin_c = sin_sin_sin + sin_cos_cos;
282 double t_sin_d = cos_sin_cos - cos_cos_sin;
283 double t_cos_a = cos_cos_cos - cos_sin_sin;
284 double t_cos_b = sin_sin_cos + sin_cos_sin;
285 double t_cos_c = cos_cos_cos + cos_sin_sin;
286 double t_cos_d = sin_sin_cos - sin_cos_sin;
287
288 struct sin_cos_angle t_angle[4] = {
289 tri_area_quarter_angle(
290 (struct sin_cos_angle){.sin = - t_sin_a + t_sin_b,
291 .cos = + t_cos_a - t_cos_b}),
292 tri_area_quarter_angle(
293 (struct sin_cos_angle){.sin = + t_sin_a + t_sin_b,
294 .cos = + t_cos_a + t_cos_b}),
295 tri_area_quarter_angle(
296 (struct sin_cos_angle){.sin = + t_sin_c - t_sin_d,
297 .cos = + t_cos_c + t_cos_d}),
298 tri_area_quarter_angle(
299 (struct sin_cos_angle){.sin = + t_sin_c + t_sin_d,
300 .cos = + t_cos_c - t_cos_d})};
301
302 double t = (t_angle[0].sin*t_angle[1].sin*t_angle[2].sin*t_angle[3].sin) /
303 (t_angle[0].cos*t_angle[1].cos*t_angle[2].cos*t_angle[3].cos);
304
305 return fabs(4.0 * atan(sqrt(fabs(t))));
306 }
307
308 /* ----------------------------------- */
309
compute_norm_vector(double a[],double b[],double norm[])310 static inline int compute_norm_vector(double a[], double b[], double norm[]) {
311
312 crossproduct_ld(a, b, norm);
313
314 double scale = sqrt(norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2]);
315
316 if (scale <= yac_angle_tol) return 1;
317
318 scale = 1.0 / scale;
319
320 norm[0] *= scale;
321 norm[1] *= scale;
322 norm[2] *= scale;
323
324 return 0;
325 }
326
XYZtoLon(double a[3])327 static double inline XYZtoLon(double a[3]) {
328 return atan2(a[1] , a[0]);
329 }
330
331 static double
lat_edge_correction(double base_vec[3],double a[3],double b[3])332 lat_edge_correction(double base_vec[3], double a[3], double b[3]) {
333
334 double const tol = 1e-8;
335
336 if (fabs(a[2] - b[2]) > tol)
337 yac_internal_abort_message("ERROR: latitude of both corners is not identical\n",
338 __FILE__, __LINE__);
339
340 double h = fabs(a[2]);
341
342 // if we are at the equator or at a pole
343 if (h < tol || fabs(1.0 - h) < tol)
344 return 0.0;
345
346 double lat_area = fabs((1.0 - h) * get_angle(XYZtoLon(a), XYZtoLon(b)));
347
348 double pole[3] = {0, 0, (a[2] >= 0.0)?1.0:-1.0};
349 double gc_area = tri_area(a, b, pole);
350
351 double correction = MAX(lat_area - gc_area, 0.0);
352
353 double middle_lat[3] = {a[0]+b[0], a[1]+b[1], a[2]};
354 double scale = sqrt(middle_lat[0]*middle_lat[0]+middle_lat[1]*middle_lat[1]);
355
356 if (fabs(scale) < 1e-18)
357 yac_internal_abort_message("internal error", __FILE__, __LINE__);
358
359 scale = sqrt(1.0 - a[2]*a[2]) / scale;
360
361 middle_lat[0] *= scale;
362 middle_lat[1] *= scale;
363
364 double norm_ab[3];
365
366 // if the angle between a and b is to small to compute a norm vector
367 if (compute_norm_vector(a, b, norm_ab)) return 0.0;
368
369 double scalar_base = scalar_product(norm_ab, base_vec);
370 double scalar_middle_lat = scalar_product(norm_ab, middle_lat);
371
372 // if the base point is on the same plane as a and b
373 if (fabs(scalar_base) < 1e-11) {
374
375 double norm_middle[3];
376
377 if (compute_norm_vector(middle_lat, pole, norm_middle)) return 0.0;
378
379 double scalar_a = scalar_product(norm_middle, a);
380
381 if (scalar_a > 0)
382 return correction;
383 else
384 return - correction;
385
386 } else {
387
388 if (scalar_middle_lat < 0)
389 return correction;
390 else
391 return - correction;
392 }
393 }
394
395
396
yac_pole_area(struct grid_cell cell)397 double yac_pole_area(struct grid_cell cell) {
398
399 size_t const M = cell.num_corners;
400
401 double coordinates_x[M];
402 double coordinates_y[M];
403
404 for (size_t i = 0; i < M; ++i)
405 XYZtoLL(cell.coordinates_xyz[i], &(coordinates_x[i]), &(coordinates_y[i]));
406
407 double area = 0.0;
408
409 if (M < 2) return 0.0;
410
411 int closer_to_south_pole = coordinates_y[0] < 0;
412
413 double pole_vec[3] = {0, 0, (closer_to_south_pole)?-1.0:1.0};
414
415 // it would also be possible to use the equator instead
416 // of the poles as the baseline
417 // this could be used as special case for cell close
418 // the equator (were the other method is probably most
419 // inaccurate)
420
421 for (size_t i = 0; i < M; ++i) {
422
423 // if one of the points it at the pole
424 if (fabs(fabs(coordinates_y[i]) - M_PI_2) < 1e-12) continue;
425 if (fabs(fabs(coordinates_y[(i+1)%M]) - M_PI_2) < 1e-12) {
426 ++i; // we can skip the next edge
427 continue;
428 }
429
430 if (cell.edge_type[i] == GREAT_CIRCLE || cell.edge_type[i] == LON_CIRCLE) {
431
432 double * a;
433 double * b;
434
435 a = cell.coordinates_xyz[i];
436 b = cell.coordinates_xyz[(i+1)%M];
437
438 double edge_direction = a[0]*b[1]-a[1]*b[0]; // 3. component of cross product
439
440 // if the edge is nearly on a circle of longitude
441 if (fabs(edge_direction) < 1e-12) continue;
442
443 double tmp_area = tri_area(a, b, pole_vec);
444
445 // or the other way round
446 if (edge_direction > 0) area -= tmp_area;
447 else area += tmp_area;
448
449 } else if (cell.edge_type[i] == LAT_CIRCLE) {
450
451 // the area of a sphere cap is:
452 // A = 2 * PI * r * h (where r == 1 and h == 1 - cos(d_lat))
453 // scaled with the longitude angle this is:
454 // A' = (d_lon / (2 * PI)) * A
455 // => A' = d_lon * (1 - cos(d_lat))
456
457 double d_lon = get_angle(coordinates_x[i], coordinates_x[(i+1)%M]);
458 double d_lat = M_PI_2;
459
460 if (closer_to_south_pole)
461 d_lat += coordinates_y[i];
462 else
463 d_lat -= coordinates_y[i];
464
465 double h = 1 - cos(d_lat);
466
467 area += d_lon * h;
468
469 } else {
470
471 yac_internal_abort_message("ERROR: unsupported edge type\n", __FILE__, __LINE__);
472 }
473 }
474 // return fabs(area) * YAC_EARTH_RADIUS * YAC_EARTH_RADIUS;
475 return fabs(area);
476 }
477
yac_planar_3dcell_area(struct grid_cell cell)478 double yac_planar_3dcell_area (struct grid_cell cell) {
479
480 /*
481 * source code is originally based on http://gaim.umbc.edu/2010/06/03/polygon-area/
482 *
483 */
484
485 double norm[3] = {0,0,0};
486 size_t M = cell.num_corners;
487
488 if (M < 3) return 0.0;
489
490 for (size_t i0 = M - 1, i1 = 0; i1 < M; i0 = i1, ++i1) {
491 norm[0] += cell.coordinates_xyz[i0][1]*cell.coordinates_xyz[i1][2] -
492 cell.coordinates_xyz[i1][1]*cell.coordinates_xyz[i0][2];
493 norm[1] += cell.coordinates_xyz[i0][2]*cell.coordinates_xyz[i1][0] -
494 cell.coordinates_xyz[i1][2]*cell.coordinates_xyz[i0][0];
495 norm[2] += cell.coordinates_xyz[i0][0]*cell.coordinates_xyz[i1][1] -
496 cell.coordinates_xyz[i1][0]*cell.coordinates_xyz[i0][1];
497 };
498
499 return 0.5 * sqrt(norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2]);
500 }
501
502 /*
503 * source code is originally based on code by Robert Oehmke of Earth System Modeling
504 * Framework (www.earthsystemmodeling.org)
505 *
506 * it has be extended to support YAC data structures and two types of
507 * grid cell edges (great circle and circle of latitude)
508 *
509 */
yac_huiliers_area(struct grid_cell cell)510 double yac_huiliers_area (struct grid_cell cell) {
511
512 size_t M = cell.num_corners;
513
514 if (M < 2) return 0.0;
515
516 int lat_flag = 0;
517
518 for (size_t i = 0; i < M; i++)
519 lat_flag |= cell.edge_type[i] == LAT_CIRCLE;
520
521 if (M == 3 && !lat_flag)
522 return fabs(tri_area(cell.coordinates_xyz[0],
523 cell.coordinates_xyz[1],
524 cell.coordinates_xyz[2]));
525
526 // sum areas around cell
527 double area = 0.0;
528
529 for (size_t i = 2; i < M; ++i) {
530
531 double tmp_area = tri_area(cell.coordinates_xyz[0],
532 cell.coordinates_xyz[i-1],
533 cell.coordinates_xyz[i]);
534
535 double norm[3];
536
537 if (compute_norm_vector(cell.coordinates_xyz[i-1],
538 cell.coordinates_xyz[i], norm)) continue;
539
540 double scalar_base = scalar_product(norm, cell.coordinates_xyz[0]);
541
542 if (scalar_base > 0) area += tmp_area;
543 else area -= tmp_area;
544 }
545
546 // if there is at least one latitude circle edge
547 if (lat_flag) {
548
549 for (size_t i = 0; i < M; ++i) {
550
551 if (cell.edge_type[i] == LAT_CIRCLE) {
552
553 size_t i_ = (i+1)%cell.num_corners;
554
555 area += lat_edge_correction(cell.coordinates_xyz[0],
556 cell.coordinates_xyz[i],
557 cell.coordinates_xyz[i_]);
558 }
559 }
560 }
561
562 return fabs(area);
563 }
564
565 /* ----------------------------------- */
566
567 // double yac_partial_area ( double a_lon, double a_lat,
568 // double b_lon, double b_lat,
569 // double c_lon, double c_lat ) {
570
571 // double theta;
572 // double angle_f;
573 // double angle_b;
574
575 // angle_f = inner_angle ( a_lat, a_lon, b_lat, b_lon );
576 // angle_b = inner_angle ( a_lat, a_lon, c_lat, c_lon );
577
578 // theta = angle_b - angle_f;
579
580 // if ( theta < 0.0 ) theta = theta + M_PI + M_PI;
581
582 // return theta;
583 // }
584
585 /* ----------------------------------- */
586
inner_angle(double plat,double plon,double qlat,double qlon)587 static inline double inner_angle ( double plat, double plon,
588 double qlat, double qlon ) {
589
590 double sin_qplon = sin(qlon-plon);
591 double cos_qlat = cos(qlat);
592 double t = sin_qplon * cos_qlat;
593
594 double sin_plat = sin(plat);
595 double sin_qlat = sin(qlat);
596 double cos_plat = cos(plat);
597 double cos_qplon = cos(qlon-plon);
598 double b = sin_qlat * cos_plat - cos_qlat * sin_plat * cos_qplon;
599
600 return atan2(b,t);
601 }
602
603 /* ----------------------------------- */
604
scalar_product(double a[],double b[])605 static inline double scalar_product(double a[], double b[]) {
606 return (a[0] * b[0] + a[1] * b[1] + a[2] * b[2]);
607 }
608
609 /* ----------------------------------- */
610
611