Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update SphericalVector to work with StructuredColumns as source functionspace. #175

Merged
merged 7 commits into from
Mar 6, 2024
Merged

Update SphericalVector to work with StructuredColumns as source functionspace. #175

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
890d5ca
Added workaround to avoid FixupHalos in interpolation test.
odlomax Mar 4, 2024
cf5c630
Realised StructuredColumns weight corrections cancel due to Euler ide…
odlomax Mar 4, 2024
ad0c800
Added great-circle course tests for points on the pole.
odlomax Mar 4, 2024
6a3f0c9
Added great-circle course special case where lonlat1 == lonlat2.
odlomax Mar 4, 2024
dc24f5a
Removed unecessary headers.
odlomax Mar 4, 2024
ffd61b3
Fixed mangled copyright.
odlomax Mar 4, 2024
27ddc3b
Fixup: atlas/util/Point.h was already included
wdeconinck Mar 6, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 2 additions & 7 deletions src/atlas/interpolation/method/sphericalvector/SphericalVector.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,13 +78,8 @@ void SphericalVector::do_setup(const FunctionSpace& source,
option::variables(2));
auto targetLonLats = target_.createField<double>(option::name("lonlat") |
option::variables(2));
sourceLonLats.array().copy(source_.lonlat());
targetLonLats.array().copy(target_.lonlat());
sourceLonLats.haloExchange();
targetLonLats.haloExchange();

const auto sourceLonLatsView = array::make_view<double, 2>(sourceLonLats);
const auto targetLonLatsView = array::make_view<double, 2>(targetLonLats);
const auto sourceLonLatsView = array::make_view<double, 2>(source_.lonlat());
const auto targetLonLatsView = array::make_view<double, 2>(target_.lonlat());

const auto unitSphere = geometry::UnitSphere{};

Expand Down
116 changes: 60 additions & 56 deletions src/atlas/util/Geometry.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
* In applying this licence, ECMWF does not waive the privileges and immGeometryies
* In applying this licence, ECMWF does not waive the privileges and immunities
Copy link
Member

@wdeconinck wdeconinck Mar 5, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for keeping our immunities intact 🙃

* granted to it by virtue of its status as an intergovernmental organisation
* nor does it submit to any jurisdiction.
*/
Expand All @@ -12,9 +12,6 @@

#include <cmath>

#include "eckit/geometry/Point2.h"
#include "eckit/geometry/Point3.h"

#include "atlas/library/config.h"
#include "atlas/runtime/Exception.h"
#include "atlas/util/Constants.h"
Expand All @@ -25,16 +22,15 @@ namespace geometry {
namespace detail {
void GeometrySphere::lonlat2xyz(const Point2& lonlat, Point3& xyz) const {
#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 24, 0)
Sphere::convertSphericalToCartesian(radius_, lonlat, xyz, 0., true);
Sphere::convertSphericalToCartesian(radius_, lonlat, xyz, 0., true);
#else
Sphere::convertSphericalToCartesian(radius_, lonlat, xyz);
Sphere::convertSphericalToCartesian(radius_, lonlat, xyz);
#endif
}
void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const {
Sphere::convertCartesianToSpherical(radius_, xyz, lonlat);
Sphere::convertCartesianToSpherical(radius_, xyz, lonlat);
}


/// @brief Calculate great-cricle course between points
///
/// @details Calculates the direction (clockwise from north) of a great-circle
Expand All @@ -45,7 +41,10 @@ void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const {
/// @ref https://en.wikipedia.org/wiki/Great-circle_navigation
///
std::pair<double, double> greatCircleCourse(const Point2& lonLat1,
const Point2& lonLat2) {
const Point2& lonLat2) {
if (lonLat1 == lonLat2) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be using bool points_equal(const Point2&, const Point2&)

I have a more consistent way to compare these in an upcoming PR (a refactoring of the geometry library) and I've already included the "course" methods there, but they will need updating to your improved version, when that time comes :-)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's correct, that it should use "points_equal". But that already happens. See https://github.com/ecmwf/atlas/blob/develop/src/atlas/util/Point.h#L32-L34.

return std::make_pair(0., 0.);
}

const auto lambda1 = lonLat1[0] * util::Constants::degreesToRadians();
const auto lambda2 = lonLat2[0] * util::Constants::degreesToRadians();
Expand All @@ -71,71 +70,76 @@ std::pair<double, double> greatCircleCourse(const Point2& lonLat1,
alpha2 * util::Constants::radiansToDegrees());
};


} // namespace detail
} // namespace geometry
} // namespace detail
} // namespace geometry

extern "C" {
// ------------------------------------------------------------------
// C wrapper interfaces to C++ routines
Geometry::Implementation* atlas__Geometry__new_name(const char* name) {
Geometry::Implementation* geometry;
{
Geometry handle{std::string{name}};
geometry = handle.get();
geometry->attach();
}
geometry->detach();
return geometry;
Geometry::Implementation* geometry;
{
Geometry handle{std::string{name}};
geometry = handle.get();
geometry->attach();
}
geometry->detach();
return geometry;
}
geometry::detail::GeometryBase* atlas__Geometry__new_radius(const double radius) {
Geometry::Implementation* geometry;
{
Geometry handle{radius};
geometry = handle.get();
geometry->attach();
}
geometry->detach();
return geometry;
geometry::detail::GeometryBase* atlas__Geometry__new_radius(
const double radius) {
Geometry::Implementation* geometry;
{
Geometry handle{radius};
geometry = handle.get();
geometry->attach();
}
geometry->detach();
return geometry;
}
void atlas__Geometry__delete(Geometry::Implementation* This) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
delete This;
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
delete This;
}
void atlas__Geometry__xyz2lonlat(Geometry::Implementation* This, const double x, const double y, const double z,
double& lon, double& lat) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
PointLonLat lonlat;
This->xyz2lonlat(PointXYZ{x, y, z}, lonlat);
lon = lonlat.lon();
lat = lonlat.lat();
void atlas__Geometry__xyz2lonlat(Geometry::Implementation* This, const double x,
const double y, const double z, double& lon,
double& lat) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
PointLonLat lonlat;
This->xyz2lonlat(PointXYZ{x, y, z}, lonlat);
lon = lonlat.lon();
lat = lonlat.lat();
}
void atlas__Geometry__lonlat2xyz(Geometry::Implementation* This, const double lon, const double lat, double& x,
void atlas__Geometry__lonlat2xyz(Geometry::Implementation* This,
const double lon, const double lat, double& x,
double& y, double& z) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
PointXYZ xyz;
This->lonlat2xyz(PointLonLat{lon, lat}, xyz);
x = xyz.x();
y = xyz.y();
z = xyz.z();
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
PointXYZ xyz;
This->lonlat2xyz(PointLonLat{lon, lat}, xyz);
x = xyz.x();
y = xyz.y();
z = xyz.z();
}
double atlas__Geometry__distance_lonlat(Geometry::Implementation* This, const double lon1, const double lat1,
double atlas__Geometry__distance_lonlat(Geometry::Implementation* This,
const double lon1, const double lat1,
const double lon2, const double lat2) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->distance(PointLonLat{lon1, lat1}, PointLonLat{lon2, lat2});
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->distance(PointLonLat{lon1, lat1}, PointLonLat{lon2, lat2});
}
double atlas__Geometry__distance_xyz(Geometry::Implementation* This, const double x1, const double y1, const double z1,
const double x2, const double y2, const double z2) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->distance(PointXYZ{x1, y1, z1}, PointXYZ{x2, y2, z2});
double atlas__Geometry__distance_xyz(Geometry::Implementation* This,
const double x1, const double y1,
const double z1, const double x2,
const double y2, const double z2) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->distance(PointXYZ{x1, y1, z1}, PointXYZ{x2, y2, z2});
}
double atlas__Geometry__radius(Geometry::Implementation* This) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->radius();
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->radius();
}
double atlas__Geometry__area(Geometry::Implementation* This) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->area();
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->area();
}
}
// ------------------------------------------------------------------
Expand Down
Loading
Loading