Commit 43343ac3 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

clipping: update to YAC version 1.0.3

parent 55476b69
2015-10-26 Uwe Schulzweida
* using CDI library version 1.6.9
* Version 1.6.9 released
* clipping: update to YAC version 1.0.3
2015-03-27 Uwe Schulzweida
* Fldstat: added parameter noweights to use constant grid cell area weights
......
......@@ -46,7 +46,7 @@
// an area of 20m x 20m on the Earth Surface is equivalent to an area on the unit sphere:
double area_tol () {
double yac_area_tol () {
return 0.02 / EarthRadius;
}
......@@ -55,13 +55,13 @@ static inline double scalar_product(double a[], double b[]);
static inline double inner_angle ( double plat, double plon,
double qlon, double qlat );
static double partial_area ( double a_lon, double a_lat,
static double yac_partial_area ( double a_lon, double a_lat,
double b_lon, double b_lat,
double c_lon, double c_lat );
/* ----------------------------------- */
double cell_approx_area ( struct grid_cell cell ) {
double yac_cell_approx_area ( struct grid_cell cell ) {
/* adopted from Robert.G. Chamberlain and William.H. Duquette */
......@@ -82,7 +82,7 @@ double cell_approx_area ( struct grid_cell cell ) {
/* ----------------------------------- */
double triangle_area ( struct grid_cell cell ) {
double yac_triangle_area ( struct grid_cell cell ) {
/* taken from the ICON code, mo_base_geometry.f90
provided by Luis Kornblueh, MPI-M. */
......@@ -165,7 +165,7 @@ double triangle_area ( struct grid_cell cell ) {
/* ----------------------------------- */
double cell_area ( struct grid_cell cell ) {
double yac_cell_area ( struct grid_cell cell ) {
/* generalised version based on the ICON code, mo_base_geometry.f90
provided by Luis Kornblueh, MPI-M. */
......@@ -246,7 +246,7 @@ double cell_area ( struct grid_cell cell ) {
/* ----------------------------------- */
double girards_area ( struct grid_cell cell ) {
double yac_girards_area ( struct grid_cell cell ) {
/* Bevis and Cambareri, 1987
......@@ -266,11 +266,11 @@ double girards_area ( struct grid_cell cell ) {
double * theta = malloc ( M * sizeof(theta[0]) );
for ( m = 0; m < M; m++ ) {
theta[m] = partial_area(cell.coordinates_x[(m+1)%M], cell.coordinates_y[(m+1)%M],
theta[m] = yac_partial_area(cell.coordinates_x[(m+1)%M], cell.coordinates_y[(m+1)%M],
cell.coordinates_x[(m+2)%M], cell.coordinates_y[(m+2)%M],
cell.coordinates_x[m%M], cell.coordinates_y[m%M]);
if ( theta[m] < tol ) {
area = cell3d_area ( cell );
area = yac_planar_3dcell_area ( cell );
return area;
}
}
......@@ -288,111 +288,13 @@ double girards_area ( struct grid_cell cell ) {
free ( theta );
if ( area <= 0.0 ) area = cell3d_area ( cell );
if ( area <= 0.0 ) area = yac_planar_3dcell_area ( cell );
return area;
}
/* ----------------------------------- */
double cell3d_area( struct grid_cell cell ) {
/* http://geomalgorithms.com/a01-_area.html */
int const M = cell.num_corners;
double area = 0.0;
if (M < 3) return 0.0; // a degenerated cell
double an, ax, ay, az; // abs value of normal and its coords
int coord; // coord to ignore: 1=x[1], 2=x[2], 3=x[3]
double * V[cell.num_corners];
double Norm[3];
double edge[2][3];
for (int i = 0; i < cell.num_corners; ++i)
V[i] = cell.coordinates_xyz + i * 3;
/* compute normal vector */
edge[0][0] = V[0][0] - V[1][0];
edge[0][1] = V[0][1] - V[1][1];
edge[0][2] = V[0][2] - V[1][2];
edge[1][0] = V[0][0] - V[2][0];
edge[1][1] = V[0][1] - V[2][1];
edge[1][2] = V[0][2] - V[2][2];
crossproduct_ld(edge[0], edge[1], Norm);
/* select largest abs coordinate to ignore for projection */
ax = (Norm[0]>0 ? Norm[0] : -Norm[0]); // abs x-coord
ay = (Norm[1]>0 ? Norm[1] : -Norm[1]); // abs y-coord
az = (Norm[2]>0 ? Norm[2] : -Norm[2]); // abs z-coord
coord = 3; // ignore z-coord x[2]
if (ax > ay) {
if (ax > az) coord = 1; // ignore x-coord x[0]
}
else if (ay > az) coord = 2; // ignore y-coord x[1]
/* compute area of the 2D projection
fabs is used to remain independent from
cyclic or anticyclic ordering of vertices */
for (int i = 1; i<M; i++) {
switch (coord) {
case 1:
area += fabs(V[i%M][1] * (V[(i+1)%M][2] - V[(i-1)%M][2]));
continue;
case 2:
area += fabs(V[i%M][0] * (V[(i+1)%M][2] - V[(i-1)%M][2]));
continue;
case 3:
area += fabs(V[i%M][0] * (V[(i+1)%M][1] - V[(i-1)%M][1]));
continue;
}
}
switch (coord) { // wrap-around term
case 1:
area += fabs(V[0][1] * (V[1][2] - V[M-1][2]));
break;
case 2:
area += fabs(V[0][0] * (V[1][2] - V[M-1][2]));
break;
case 3:
area += fabs(V[0][0] * (V[1][1] - V[M-1][1]));
break;
}
/* scale to get area before projection */
an = sqrt( ax*ax + ay*ay + az*az); // length of normal vector
switch (coord) {
case 1:
if ( ax > 0 )
area *= (an / (2*ax));
else
area = 0.0;
break;
case 2:
if ( ay > 0 )
area *= (an / (2*ay));
else
area = 0.0;
break;
case 3:
if ( az > 0 )
area *= (an / (2*az));
else
area = 0.0;
break;
}
// return area * EarthRadius * EarthRadius;
return area;
}
/** area of a spherical triangle based on L'Huilier's Theorem
*
* source code is taken from code by Robert Oehmke of Earth System Modeling
......@@ -454,6 +356,8 @@ tri_area(double u[3], double v[3], double w[3]) {
return fabs(4.0 * atan(sqrt(fabs(t))));
}
/* ----------------------------------- */
static inline int compute_norm_vector(double a[], double b[], double norm[]) {
crossproduct_ld(a, b, norm);
......@@ -478,8 +382,8 @@ lat_edge_correction(double base_point[3], double a[3], double b[3],
double const tol = 1e-8;
if (fabs(a[2] - b[2]) > tol)
abort_message("ERROR: latitude of both corners is not identical\n",
__FILE__, __LINE__);
yac_internal_abort_message("ERROR: latitude of both corners is not identical\n",
__FILE__, __LINE__);
double h = fabs(a[2]);
......@@ -497,7 +401,7 @@ lat_edge_correction(double base_point[3], double a[3], double b[3],
double middle_lat[3] = {a[0]+b[0], a[1]+b[1], a[2]};
double scale = sqrt(middle_lat[0]*middle_lat[0]+middle_lat[1]*middle_lat[1]);
if (scale == 0) abort_message("internal error", __FILE__, __LINE__);
if (scale == 0) yac_internal_abort_message("internal error", __FILE__, __LINE__);
scale = sqrt(1.0 - a[2]*a[2]) / scale;
......@@ -536,7 +440,7 @@ lat_edge_correction(double base_point[3], double a[3], double b[3],
}
}
double pole_area ( struct grid_cell cell ) {
double yac_pole_area ( struct grid_cell cell ) {
double area = 0.0;
......@@ -604,13 +508,35 @@ double pole_area ( struct grid_cell cell ) {
} else {
abort_message("ERROR: unsupported edge type\n", __FILE__, __LINE__);
yac_internal_abort_message("ERROR: unsupported edge type\n", __FILE__, __LINE__);
}
}
// return fabs(area) * EarthRadius * EarthRadius;
return fabs(area);
}
double yac_planar_3dcell_area (struct grid_cell cell) {
/*
* source code is originally based on http://gaim.umbc.edu/2010/06/03/polygon-area/
*
*/
double area = 0.0;
double norm[3] = {0,0,0};
if (cell.num_corners < 3) return area;
for ( int i0=cell.num_corners-1, i1=0; i1<cell.num_corners; i0=i1, ++i1) {
norm[0] += cell.coordinates_xyz[1+i0*3]*cell.coordinates_xyz[2+i1*3] - cell.coordinates_xyz[1+i1*3]*cell.coordinates_xyz[2+i0*3];
norm[1] += cell.coordinates_xyz[2+i0*3]*cell.coordinates_xyz[0+i1*3] - cell.coordinates_xyz[2+i1*3]*cell.coordinates_xyz[0+i0*3];
norm[2] += cell.coordinates_xyz[0+i0*3]*cell.coordinates_xyz[1+i1*3] - cell.coordinates_xyz[0+i1*3]*cell.coordinates_xyz[1+i0*3];
};
return area = 0.5 * sqrt(norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2] );
}
/*
* source code is originally based on code by Robert Oehmke of Earth System Modeling
* Framework (www.earthsystemmodeling.org)
......@@ -619,7 +545,7 @@ double pole_area ( struct grid_cell cell ) {
* grid cell edges (great circle and circle of latitude)
*
*/
double huiliers_area(struct grid_cell cell) {
double yac_huiliers_area (struct grid_cell cell) {
if (cell.num_corners < 2) return 0;
......@@ -680,9 +606,9 @@ double huiliers_area(struct grid_cell cell) {
/* ----------------------------------- */
double partial_area ( double a_lon, double a_lat,
double b_lon, double b_lat,
double c_lon, double c_lat ) {
double yac_partial_area ( double a_lon, double a_lat,
double b_lon, double b_lat,
double c_lon, double c_lat ) {
double theta;
double angle_f;
......
......@@ -42,7 +42,7 @@
*
**/
double area_tol();
double yac_area_tol();
/** \brief Calculate the area of a cell on a unit sphere
*
......@@ -96,7 +96,7 @@ double area_tol();
*
**/
double cell_approx_area ( struct grid_cell cell );
double yac_cell_approx_area ( struct grid_cell cell );
/** \brief Calculate the area of a triangle on a unit sphere
*
......@@ -130,7 +130,7 @@ double cell_approx_area ( struct grid_cell cell );
*
**/
double triangle_area ( struct grid_cell cell );
double yac_triangle_area ( struct grid_cell cell );
/** \brief Calculate the area of a cell on a unit sphere
*
......@@ -159,7 +159,7 @@ double triangle_area ( struct grid_cell cell );
*
**/
double cell_area ( struct grid_cell cell );
double yac_cell_area ( struct grid_cell cell );
/** \brief Calculate the area of a cell on a unit sphere
*
......@@ -183,12 +183,9 @@ double cell_area ( struct grid_cell cell );
*
* S = [ Theta - (n-2) * pi ] * R*R
*
* \todo
* - implement an alternative for very small elements
*
* This algorithm does not work for very small elements. In this case negative
* areas are delivered which are set to zero currently. An alternative would be
* to transform the coordinates into the cartesian space, ignore the curvature
* to transform the coordinates into the s space, ignore the curvature
* of the edges (as they can be neglected for small elements) and calculate the
* area of a plane triangle.
*
......@@ -201,7 +198,7 @@ double cell_area ( struct grid_cell cell );
* @return area of the cell
*
**/
double girards_area ( struct grid_cell cell );
double yac_girards_area ( struct grid_cell cell );
/** \brief Calculate the area of a cell in a 3d plane on a unit sphere
*
......@@ -219,24 +216,21 @@ double girards_area ( struct grid_cell cell );
*
**/
double cell3d_area( struct grid_cell cell );
double yac_pole_area ( struct grid_cell cell );
/** \brief Calculate the area of a cell on a unit sphere
*
* based on:
/**
* \brief Area calculation on a unit sphere of a planar polygon in 3D
*
* Robert.G. Chamberlain, William.H. Duquette
* Jet Propulsion Laboratory
* Some Algorithms for Polygons on a Sphere.
* Association of American Geographers Annual Meeting San Francisco, California
* 17 - 21 April 2007
* http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/40409/1/07-03.pdf
* (http://gaim.umbc.edu/2010/06/03/polygon-area)\n
*
* uses different formula to compute the area of the pole triangles
* This area calculation works for any planar polygon (concave or convex)
* with non-intersecting edges in 3D. It requires vertex coordinates in
* Carthesian space. In our case this is applicable for very
* small elements on the sphere.
*
* \todo enable routine to be able to handle latitude circle edges
*/
double pole_area ( struct grid_cell cell );
double yac_planar_3dcell_area (struct grid_cell cell);
/**
* \brief Area calculation on a unit sphere taken from ESMF based on L'Huilier's Theorem
......@@ -248,8 +242,10 @@ double pole_area ( struct grid_cell cell );
* The cell in split up into triangles that all have one corner in common,
* then the area for each of the triangle is computed and summed up to build
* the area of the cell. L'Huilier's Theorem is used to compute the area of
* the triangles.
* the triangles. This seems to be sufficiently accurate for elements on the
* Earth surface with edge lengths of approx. 100 m and gives results comparable
* to our implementation of Huilier's algorithm for edge lengths of up to 1 km.
*/
double huiliers_area(struct grid_cell cell);
double yac_huiliers_area(struct grid_cell cell);
#endif // AREA_H
......@@ -55,7 +55,7 @@ enum cell_type {
struct point_list_element {
double vec_coords[3];
enum edge_type edge_type; // type of edge with next corner
enum yac_edge_type edge_type; // type of edge with next corner
int to_be_removed;
struct point_list_element * next;
};
......@@ -92,10 +92,10 @@ static enum cell_type get_cell_type(struct grid_cell target_cell);
/* ------------------------- */
void compute_overlap_areas(unsigned N,
struct grid_cell * source_cell,
struct grid_cell target_cell,
double * partial_areas) {
void yac_compute_overlap_areas (unsigned N,
struct grid_cell * source_cell,
struct grid_cell target_cell,
double * partial_areas) {
static struct grid_cell * overlap_buffer = NULL;
static unsigned overlap_buffer_size = 0;
......@@ -110,17 +110,17 @@ void compute_overlap_areas(unsigned N,
for (; old_overlap_buffer_size < overlap_buffer_size;
++old_overlap_buffer_size)
init_grid_cell(overlap_buffer + old_overlap_buffer_size);
yac_init_grid_cell(overlap_buffer + old_overlap_buffer_size);
}
/* Do the clipping and get the cell for the overlapping area */
cell_clipping ( N, source_cell, target_cell, overlap_buffer);
yac_cell_clipping ( N, source_cell, target_cell, overlap_buffer);
/* Get the partial areas for the overlapping regions */
for (unsigned n = 0; n < N; n++) {
partial_areas[n] = huiliers_area (overlap_buffer[n]);
partial_areas[n] = yac_huiliers_area (overlap_buffer[n]);
// we cannot use pole_area because it is rather inaccurate for great circle
// edges that are nearly circles of longitude
//partial_areas[n] = pole_area (overlap_buffer[n]);
......@@ -134,32 +134,32 @@ void compute_overlap_areas(unsigned N,
/* ------------------------- */
void compute_concave_overlap_areas (unsigned N,
struct grid_cell * source_cell,
struct grid_cell target_cell,
double * target_node_x,
double * target_node_y,
double * partial_areas) {
void yac_compute_concave_overlap_areas (unsigned N,
struct grid_cell * source_cell,
struct grid_cell target_cell,
double * target_node_x,
double * target_node_y,
double * partial_areas) {
enum cell_type target_cell_type;
if ( target_cell.num_corners > 3 )
target_cell_type = get_cell_type (target_cell);
if ( target_cell.num_corners < 4 || target_cell_type == LON_LAT_CELL ) {
compute_overlap_areas (N, source_cell, target_cell, partial_areas);
yac_compute_overlap_areas (N, source_cell, target_cell, partial_areas);
return;
}
if ( target_node_x == NULL || target_node_y == NULL )
abort_message("ERROR: missing target point coordinates "
"(x_coordinates == NULL || y_coordinates == NULL)",
__FILE__, __LINE__);
yac_internal_abort_message("ERROR: missing target point coordinates "
"(x_coordinates == NULL || y_coordinates == NULL)",
__FILE__ , __LINE__);
struct grid_cell target_partial_cell =
{.coordinates_x = (double[3]){-1},
.coordinates_y = (double[3]){-1},
.coordinates_xyz = (double[3*3]){-1},
.edge_type = (enum edge_type[3]) {GREAT_CIRCLE},
.edge_type = (enum yac_edge_type[3]) {GREAT_CIRCLE},
.num_corners = 3};
static struct grid_cell * overlap_buffer = NULL;
......@@ -175,7 +175,7 @@ void compute_concave_overlap_areas (unsigned N,
for (; old_overlap_buffer_size < overlap_buffer_size;
++old_overlap_buffer_size)
init_grid_cell(overlap_buffer + old_overlap_buffer_size);
yac_init_grid_cell(overlap_buffer + old_overlap_buffer_size);
}
/* Do the clipping and get the cell for the overlapping area */
......@@ -232,12 +232,12 @@ void compute_concave_overlap_areas (unsigned N,
target_partial_cell.coordinates_xyz[1+3*2] = target_cell.coordinates_xyz[1+3*corner_b];
target_partial_cell.coordinates_xyz[2+3*2] = target_cell.coordinates_xyz[2+3*corner_b];
cell_clipping ( N, source_cell, target_partial_cell, overlap_buffer);
yac_cell_clipping ( N, source_cell, target_partial_cell, overlap_buffer);
/* Get the partial areas for the overlapping regions as sum over the partial target cells. */
for (unsigned n = 0; n < N; n++) {
partial_areas[n] += huiliers_area (overlap_buffer[n]);
partial_areas[n] += yac_huiliers_area (overlap_buffer[n]);
// we cannot use pole_area because it is rather inaccurate for great circle
// edges that are nearly circles of longitude
//partial_areas[n] = pole_area (overlap_buffer[n]);
......@@ -264,8 +264,8 @@ static void compute_norm_vector(double a[], double b[], double norm[]) {
if ((fabs(norm[0]) < tol) &&
(fabs(norm[1]) < tol) &&
(fabs(norm[2]) < tol))
abort_message("ERROR: a and b are identical -> no norm vector\n",
__FILE__, __LINE__);
yac_internal_abort_message("ERROR: a and b are identical -> no norm vector\n",
__FILE__, __LINE__);
double scale = 1.0 / sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
......@@ -310,7 +310,7 @@ static unsigned is_inside_gc(double point[], double norm_vec[]) {
dot = dotproduct(point, norm_vec);
// if the point is on the line
if (fabs(dot) <= angle_tol * 1e-3)
if (fabs(dot) <= yac_angle_tol * 1e-3)
return 2;
return dot < 0;
......@@ -320,12 +320,12 @@ static unsigned is_inside_latc(double point[], double z) {
double temp = fabs(point[2] + z);
if (fabs(1.0 - temp) <= angle_tol * 1e-3) return 2;
if (fabs(1.0 - temp) <= yac_angle_tol * 1e-3) return 2;
else return temp < 1.0;
}
static unsigned is_inside(double point[], double help_vec[],
enum edge_type edge_type,
enum yac_edge_type edge_type,
unsigned cell_points_ordering) {
unsigned ret_val = 0;
......@@ -340,7 +340,7 @@ static unsigned is_inside(double point[], double help_vec[],
ret_val = is_inside_latc(point, help_vec[2]);
break;
default:
abort_message("invalid edge type\n", __FILE__, __LINE__);
yac_internal_abort_message("invalid edge type\n", __FILE__, __LINE__);
};
if (ret_val == 2) return 2;
......@@ -425,8 +425,8 @@ static void get_edge_middle_point_gc(double a[3], double b[3],
middle[2] *= scale;
}
static void get_edge_middle_point(double a[3], double b[3],
enum edge_type edge_type, double middle[3]) {
static void get_edge_middle_point (double a[3], double b[3],
enum yac_edge_type edge_type, double middle[3]) {
switch (edge_type) {
......@@ -442,7 +442,7 @@ static void get_edge_middle_point(double a[3], double b[3],
break;
default:
abort_message("ERROR: invalid edge type\n", __FILE__, __LINE__);
yac_internal_abort_message("ERROR: invalid edge type\n", __FILE__, __LINE__);
middle[0] = -1; // program should never reach this point
middle[1] = -1;
middle[2] = -1;
......@@ -450,11 +450,11 @@ static void get_edge_middle_point(double a[3], double b[3],
}
/**
* cell clipping using SutherlandHodgman algorithm;
* cell clipping using Sutherland-Hodgman algorithm;
*/
void point_list_clipping(struct point_list * source_list, int source_ordering,
struct point_list target_list, int target_ordering,
unsigned nct, double * tgt_edge_norm_vec) {
static void point_list_clipping (struct point_list * source_list, int source_ordering,
struct point_list target_list, int target_ordering,
unsigned nct, double * tgt_edge_norm_vec) {
struct {
double * edge_norm_vec;
......@@ -516,13 +516,13 @@ void point_list_clipping(struct point_list * source_list, int source_ordering,
(prev_is_inside + curr_is_inside < 4))) {
// get intersection points
intersect = intersect_vec(prev_src_point->edge_type,
prev_src_point->vec_coords,
curr_src_point->vec_coords,
tgt_points[i].point->edge_type,
tgt_points[i].point->vec_coords,
tgt_points[i].point->next->vec_coords,
p, q);
intersect = yac_intersect_vec(prev_src_point->edge_type,
prev_src_point->vec_coords,
curr_src_point->vec_coords,
tgt_points[i].point->edge_type,
tgt_points[i].point->vec_coords,
tgt_points[i].point->next->vec_coords,
p, q);
// if both edges are on an identical great circle
if ((intersect != -1) && (intersect & (1 << 4))) {
......@@ -567,7 +567,7 @@ void point_list_clipping(struct point_list * source_list, int source_ordering,
prev_src_point:curr_src_point;
break;
default:
abort_message("invalid edge type\n", __FILE__, __LINE__);
yac_internal_abort_message("invalid edge type\n", __FILE__, __LINE__);
return;
};
......@@ -582,14 +582,14 @@ void point_list_clipping(struct point_list * source_list, int source_ordering,
if (!((prev_src_point->edge_type == LAT_CIRCLE) ^
(tgt_points[i].point->edge_type == LAT_CIRCLE))) {
abort_message("ERROR: ...this should not have happened...\n",
__FILE__, __LINE__);
yac_internal_abort_message("ERROR: ...this should not have happened...\n",
__FILE__, __LINE__);
}
if (((get_vector_angle(prev_src_point->vec_coords, p) < angle_tol) &&
(get_vector_angle(curr_src_point->vec_coords, q) < angle_tol)) ||
((get_vector_angle(prev_src_point->vec_coords, q) < angle_tol) &&
(get_vector_angle(curr_src_point->vec_coords, p) < angle_tol))) {
if (((get_vector_angle(prev_src_point->vec_coords, p) < yac_angle_tol) &&
(get_vector_angle(curr_src_point->vec_coords, q) < yac_angle_tol)) ||
((get_vector_angle(prev_src_point->vec_coords, q) < yac_angle_tol) &&
(get_vector_angle(curr_src_point->vec_coords, p) < yac_angle_tol))) {
prev_is_inside = 2;
curr_is_inside = 2;
......@@ -619,7 +619,7 @@ void point_list_clipping(struct point_list * source_list, int source_ordering,
tgt_points[i].edge_norm_vec[2]));
break;