Add a jacobian method to atlas::Projection object

src/atlas/projection/Projection.cc

@@ -50,6 +50,10 @@ void atlas::Projection::lonlat2xy( Point2& point ) const {
     return get()->lonlat2xy( point );
 }
 
+atlas::Projection::Jacobian atlas::Projection::jacobian( const PointLonLat& p ) const {
+    return get()->jacobian( p );
+}
+
 PointLonLat atlas::Projection::lonlat( const PointXY& xy ) const {
     return get()->lonlat( xy );
 }

src/atlas/projection/Projection.h

@@ -14,6 +14,7 @@
 
 #include "atlas/domain/Domain.h"
 #include "atlas/library/config.h"
+#include "atlas/projection/detail/ProjectionImpl.h"
 #include "atlas/util/ObjectHandle.h"
 
 //---------------------------------------------------------------------------------------------------------------------
@@ -47,7 +48,8 @@ class ProjectionImpl;
 
 class Projection : DOXYGEN_HIDE( public util::ObjectHandle<projection::detail::ProjectionImpl> ) {
 public:
-    using Spec = util::Config;
+    using Spec     = util::Config;
+    using Jacobian = projection::detail::ProjectionImpl::Jacobian;
 
 public:
     using Handle::Handle;
@@ -62,6 +64,8 @@ class Projection : DOXYGEN_HIDE( public util::ObjectHandle<projection::detail::P
     void lonlat2xy( double crd[] ) const;
     void lonlat2xy( Point2& ) const;
 
+    Jacobian jacobian( const PointLonLat& ) const;
+
     PointLonLat lonlat( const PointXY& ) const;
     PointXY xy( const PointLonLat& ) const;
 

src/atlas/projection/detail/LambertAzimuthalEqualAreaProjection.cc

@@ -86,6 +86,10 @@ void LambertAzimuthalEqualAreaProjection::xy2lonlat( double crd[] ) const {
 }
 
 
+ProjectionImpl::Jacobian LambertAzimuthalEqualAreaProjection::jacobian( const PointLonLat& ) const {
+    throw_NotImplemented( "LambertAzimuthalEqualAreaProjection::jacobian", Here() );
+}
+
 LambertAzimuthalEqualAreaProjection::Spec LambertAzimuthalEqualAreaProjection::spec() const {
     Spec proj;
     proj.set( "type", static_type() );

src/atlas/projection/detail/LambertAzimuthalEqualAreaProjection.h

@@ -30,6 +30,8 @@ class LambertAzimuthalEqualAreaProjection final : public ProjectionImpl {
     void xy2lonlat( double crd[] ) const override;
     void lonlat2xy( double crd[] ) const override;
 
+    Jacobian jacobian( const PointLonLat& ) const override;
+
     bool strictlyRegional() const override { return true; }
     RectangularLonLatDomain lonlatBoundingBox( const Domain& domain ) const override {
         return ProjectionImpl::lonlatBoundingBox( domain );

src/atlas/projection/detail/LambertConformalConicProjection.cc

@@ -100,6 +100,38 @@ void LambertConformalConicProjection::xy2lonlat( double crd[] ) const {
             : util::Constants::radiansToDegrees() * 2. * std::atan( std::pow( radius_ * F_ / rho, inv_n_ ) ) - 90.;
 }
 
+ProjectionImpl::Jacobian LambertConformalConicProjection::jacobian( const PointLonLat& lonlat ) const {
+    ProjectionImpl::Jacobian jac;
+
+    const double deg2rad = util::Constants::degreesToRadians();
+
+    double lat = lonlat.lat(), lon = lonlat.lon();
+    double tanlat = tan_d( lat ), coslat = cos_d( 45. + lat * 0.5 );
+
+    double rho   = radius_ * F_ * std::pow( tanlat, -n_ );
+    double theta = n_ * normalise( lon - lon0_, -180, 360 );
+
+    double costheta = cos_d( theta ), sintheta = sin_d( theta );
+
+    double x = rho * sintheta;
+    double y = rho0_ - rho * costheta;
+
+    auto cpj = [&]( const double dlon, const double dlat, double& dx, double& dy ) {
+        double drho = deg2rad * radius_ * F_ * ( -n_ ) * std::pow( tanlat, -n_ - 1 ) / ( coslat * coslat ) * dlat * 0.5;
+        double dtheta = +n_ * dlon;
+        dx            = +drho * sintheta + rho * costheta * deg2rad * dtheta;
+        dy            = -drho * costheta + rho * sintheta * deg2rad * dtheta;
+    };
+
+    const double dlat = 1.0f;
+    const double dlon = 1.0f;
+
+    cpj( dlon, 0, jac[0][0], jac[1][0] );
+    cpj( 0, dlat, jac[0][1], jac[1][1] );
+
+    return jac;
+}
+
 LambertConformalConicProjection::Spec LambertConformalConicProjection::spec() const {
     Spec spec;
     spec.set( "type", static_type() );

src/atlas/projection/detail/LambertConformalConicProjection.h

@@ -30,6 +30,8 @@ class LambertConformalConicProjection final : public ProjectionImpl {
     void xy2lonlat( double crd[] ) const override;
     void lonlat2xy( double crd[] ) const override;
 
+    Jacobian jacobian( const PointLonLat& ) const override;
+
     bool strictlyRegional() const override { return true; }
 
     RectangularLonLatDomain lonlatBoundingBox( const Domain& domain ) const override {

src/atlas/projection/detail/LonLatProjection.cc

@@ -35,6 +35,19 @@ void LonLatProjectionT<NotRotated>::xy2lonlat( double[] ) const {}
 template <>
 void LonLatProjectionT<NotRotated>::lonlat2xy( double[] ) const {}
 
+template <>
+ProjectionImpl::Jacobian LonLatProjectionT<NotRotated>::jacobian( const PointLonLat& ) const {
+    Jacobian jac;
+    jac[0] = {1.0, 0.0};
+    jac[1] = {0.0, 1.0};
+    return jac;
+}
+
+template <typename Rotation>
+ProjectionImpl::Jacobian LonLatProjectionT<Rotation>::jacobian( const PointLonLat& ) const {
+    throw_NotImplemented( "LonLatProjectionT::jacobian", Here() );
+}
+
 template <>
 RectangularLonLatDomain LonLatProjectionT<NotRotated>::lonlatBoundingBox( const Domain& domain ) const {
     return domain;

src/atlas/projection/detail/LonLatProjection.h

@@ -41,6 +41,8 @@ class LonLatProjectionT final : public ProjectionImpl {
     void xy2lonlat( double crd[] ) const override { rotation_.rotate( crd ); }
     void lonlat2xy( double crd[] ) const override { rotation_.unrotate( crd ); }
 
+    Jacobian jacobian( const PointLonLat& ) const override;
+
     bool strictlyRegional() const override { return false; }
     RectangularLonLatDomain lonlatBoundingBox( const Domain& ) const override;
 

src/atlas/projection/detail/MercatorProjection.cc

@@ -153,6 +153,11 @@ void MercatorProjectionT<Rotation>::xy2lonlat( double crd[] ) const {
     normalise_( crd );
 }
 
+template <typename Rotation>
+ProjectionImpl::Jacobian MercatorProjectionT<Rotation>::jacobian( const PointLonLat& ) const {
+    throw_NotImplemented( "MercatorProjectionT::jacobian", Here() );
+}
+
 // specification
 template <typename Rotation>
 typename MercatorProjectionT<Rotation>::Spec MercatorProjectionT<Rotation>::spec() const {

src/atlas/projection/detail/MercatorProjection.h

@@ -32,6 +32,8 @@ class MercatorProjectionT final : public ProjectionImpl {
     void xy2lonlat( double crd[] ) const override;
     void lonlat2xy( double crd[] ) const override;
 
+    Jacobian jacobian( const PointLonLat& ) const override;
+
     bool strictlyRegional() const override { return true; }  // Mercator projection cannot be used for global grids
     RectangularLonLatDomain lonlatBoundingBox( const Domain& domain ) const override {
         return ProjectionImpl::lonlatBoundingBox( domain );

src/atlas/projection/detail/ProjProjection.cc

@@ -144,6 +144,10 @@ void ProjProjection::lonlat2xy( double crd[] ) const {
 }
 
 
+ProjectionImpl::Jacobian ProjProjection::jacobian( const PointLonLat& ) const {
+    throw_NotImplemented( "ProjProjection::jacobian", Here() );
+}
+
 PointXYZ ProjProjection::xyz( const PointLonLat& lonlat ) const {
     PJ_COORD P = proj_coord( lonlat.lon(), lonlat.lat(), 0, 0 );
     P          = proj_trans( sourceToGeocentric_, PJ_FWD, P );

src/atlas/projection/detail/ProjProjection.h

@@ -55,6 +55,9 @@ class ProjProjection final : public ProjectionImpl {
 
     void xy2lonlat( double[] ) const override;
     void lonlat2xy( double[] ) const override;
+
+    Jacobian jacobian( const PointLonLat& ) const override;
+
     PointXYZ xyz( const PointLonLat& ) const override;
 
     bool strictlyRegional() const override { return false; }

src/atlas/projection/detail/ProjectionImpl.h

@@ -39,6 +39,114 @@ namespace detail {
 class ProjectionImpl : public util::Object {
 public:
     using Spec = atlas::util::Config;
+    class Jacobian : public std::array<std::array<double, 2>, 2> {
+    public:
+        static Jacobian identity() {
+            Jacobian id;
+            id[0] = {1.0, 0.0};
+            id[1] = {0.0, 1.0};
+            return id;
+        }
+
+        Jacobian inverse() const {
+            const Jacobian& jac = *this;
+            Jacobian inv;
+            double det = jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
+            inv[0][0]  = +jac[1][1] / det;
+            inv[0][1]  = -jac[0][1] / det;
+            inv[1][0]  = -jac[1][0] / det;
+            inv[1][1]  = +jac[0][0] / det;
+            return inv;
+        };
+
+        Jacobian transpose() const {
+            Jacobian tra = *this;
+            std::swap( tra[0][1], tra[1][0] );
+            return tra;
+        };
+
+        Jacobian operator-( const Jacobian& jac2 ) const {
+            const Jacobian& jac1 = *this;
+            Jacobian jac;
+            jac[0][0] = jac1[0][0] - jac2[0][0];
+            jac[0][1] = jac1[0][1] - jac2[0][1];
+            jac[1][0] = jac1[1][0] - jac2[1][0];
+            jac[1][1] = jac1[1][1] - jac2[1][1];
+            return jac;
+        }
+
+        Jacobian operator+( const Jacobian& jac2 ) const {
+            const Jacobian& jac1 = *this;
+            Jacobian jac;
+            jac[0][0] = jac1[0][0] + jac2[0][0];
+            jac[0][1] = jac1[0][1] + jac2[0][1];
+            jac[1][0] = jac1[1][0] + jac2[1][0];
+            jac[1][1] = jac1[1][1] + jac2[1][1];
+            return jac;
+        }
+
+        double norm() const {
+            const Jacobian& jac = *this;
+            return sqrt( jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1] + jac[1][0] * jac[1][0] +
+                         jac[1][1] * jac[1][1] );
+        }
+
+        Jacobian operator*( const Jacobian& jac2 ) const {
+            const Jacobian& jac1 = *this;
+            Jacobian jac;
+            jac[0][0] = jac1[0][0] * jac2[0][0] + jac1[0][1] * jac2[1][0];
+            jac[0][1] = jac1[0][0] * jac2[0][1] + jac1[0][1] * jac2[1][1];
+            jac[1][0] = jac1[1][0] * jac2[0][0] + jac1[1][1] * jac2[1][0];
+            jac[1][1] = jac1[1][0] * jac2[0][1] + jac1[1][1] * jac2[1][1];
+            return jac;
+        }
+
+        double dx_dlon() const {
+            const Jacobian& jac = *this;
+            return jac[JDX][JDLON];
+        }
+        double dy_dlon() const {
+            const Jacobian& jac = *this;
+            return jac[JDY][JDLON];
+        }
+        double dx_dlat() const {
+            const Jacobian& jac = *this;
+            return jac[JDX][JDLAT];
+        }
+        double dy_dlat() const {
+            const Jacobian& jac = *this;
+            return jac[JDY][JDLAT];
+        }
+
+        double dlon_dx() const {
+            const Jacobian& jac = *this;
+            return jac[JDLON][JDX];
+        }
+        double dlon_dy() const {
+            const Jacobian& jac = *this;
+            return jac[JDLON][JDY];
+        }
+        double dlat_dx() const {
+            const Jacobian& jac = *this;
+            return jac[JDLAT][JDX];
+        }
+        double dlat_dy() const {
+            const Jacobian& jac = *this;
+            return jac[JDLAT][JDY];
+        }
+
+    private:
+        enum
+        {
+            JDX = 0,
+            JDY = 1
+        };
+        enum
+        {
+            JDLON = 0,
+            JDLAT = 1
+        };
+    };
 
 public:
     static const ProjectionImpl* create( const eckit::Parametrisation& p );
@@ -52,6 +160,8 @@ class ProjectionImpl : public util::Object {
     virtual void xy2lonlat( double crd[] ) const = 0;
     virtual void lonlat2xy( double crd[] ) const = 0;
 
+    virtual Jacobian jacobian( const PointLonLat& ) const = 0;
+
     void xy2lonlat( Point2& ) const;
     void lonlat2xy( Point2& ) const;
 

src/atlas/projection/detail/SchmidtProjection.cc

@@ -18,6 +18,7 @@
 #include "atlas/runtime/Exception.h"
 #include "atlas/util/Config.h"
 #include "atlas/util/Constants.h"
+#include "atlas/util/UnitSphere.h"
 
 namespace {
 static double D2R( const double x ) {
@@ -65,6 +66,53 @@ void SchmidtProjectionT<Rotation>::lonlat2xy( double crd[] ) const {
         R2D( std::asin( std::cos( 2. * std::atan( c_ * std::tan( std::acos( std::sin( D2R( crd[1] ) ) ) * 0.5 ) ) ) ) );
 }
 
+template <>
+ProjectionImpl::Jacobian SchmidtProjectionT<NotRotated>::jacobian( const PointLonLat& ) const {
+    throw_NotImplemented( "SchmidtProjectionT<NotRotated>::jacobian", Here() );
+}
+
+template <typename Rotation>
+ProjectionImpl::Jacobian SchmidtProjectionT<Rotation>::jacobian( const PointLonLat& lonlat ) const {
+    double xy[2] = {lonlat.lon(), lonlat.lat()};
+
+    lonlat2xy( xy );
+
+    PointXYZ xyz, north1, north0( 0.0, 0.0, 1.0 );
+    atlas::util::UnitSphere::convertSphericalToCartesian( rotation_.northPole(), north1 );
+    atlas::util::UnitSphere::convertSphericalToCartesian( lonlat, xyz );
+
+    north1 = PointXYZ::normalize( north1 );
+
+    // Base vectors in unrotated frame
+    auto u0 = PointXYZ::normalize( PointXYZ::cross( north0, xyz ) );
+    auto v0 = PointXYZ::normalize( PointXYZ::cross( xyz, u0 ) );
+
+    // Base vectors in rotated frame
+    auto u1 = PointXYZ::normalize( PointXYZ::cross( north1, xyz ) );
+    auto v1 = PointXYZ::normalize( PointXYZ::cross( xyz, u1 ) );
+
+    double u0u1 = PointXYZ::dot( u0, u1 );
+    double v0u1 = PointXYZ::dot( v0, u1 );
+    double u0v1 = PointXYZ::dot( u0, v1 );
+    double v0v1 = PointXYZ::dot( v0, v1 );
+
+    Jacobian jac;
+
+    double zomc2 = 1.0 - 1.0 / ( c_ * c_ );
+    double zopc2 = 1.0 + 1.0 / ( c_ * c_ );
+
+    double zcosy = cos( D2R( xy[1] ) ), zsiny = sin( D2R( xy[1] ) );
+    double zcoslat = cos( D2R( lonlat.lat() ) );
+
+    double zfactor = sqrt( ( zopc2 + zsiny * zomc2 ) * ( zopc2 + zsiny * zomc2 ) / ( zopc2 * zopc2 - zomc2 * zomc2 ) );
+
+
+    jac[0] = {zcoslat * u0u1 * zfactor / zcosy, v0u1 * zfactor / zcosy};
+    jac[1] = {zcoslat * u0v1 * zfactor, v0v1 * zfactor};
+
+    return jac;
+}
+
 // specification
 template <typename Rotation>
 typename SchmidtProjectionT<Rotation>::Spec SchmidtProjectionT<Rotation>::spec() const {

src/atlas/projection/detail/SchmidtProjection.h

@@ -32,6 +32,8 @@ class SchmidtProjectionT final : public ProjectionImpl {
     void xy2lonlat( double crd[] ) const override;
     void lonlat2xy( double crd[] ) const override;
 
+    Jacobian jacobian( const PointLonLat& ) const override;
+
     bool strictlyRegional() const override { return false; }  // schmidt is global grid
     RectangularLonLatDomain lonlatBoundingBox( const Domain& domain ) const override {
         return ProjectionImpl::lonlatBoundingBox( domain );

src/tests/projection/CMakeLists.txt

@@ -9,6 +9,7 @@
 foreach(test
           test_bounding_box
           test_projection_LAEA
+          test_jacobian
           test_rotation )
 
     ecbuild_add_test( TARGET atlas_${test} SOURCES ${test}.cc LIBS atlas ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} )