diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index c558f11aa..4731266c1 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -79,6 +79,7 @@ if(ENABLE_ECL_INPUT)
     src/opm/parser/eclipse/EclipseState/Grid/Fault.cpp
     src/opm/parser/eclipse/EclipseState/Grid/FaultFace.cpp
     src/opm/parser/eclipse/EclipseState/Grid/GridDims.cpp
+    src/opm/parser/eclipse/EclipseState/Grid/MapAxes.cpp
     src/opm/parser/eclipse/EclipseState/Grid/MULTREGTScanner.cpp
     src/opm/parser/eclipse/EclipseState/Grid/NNC.cpp
     src/opm/parser/eclipse/EclipseState/Grid/Operate.cpp
@@ -632,6 +633,7 @@ if(ENABLE_ECL_INPUT)
        opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
        opm/parser/eclipse/EclipseState/Grid/BoxManager.hpp
        opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp
+       opm/parser/eclipse/EclipseState/Grid/MapAxes.hpp
        opm/parser/eclipse/EclipseState/Grid/MinpvMode.hpp
        opm/parser/eclipse/EclipseState/EndpointScaling.hpp
        opm/parser/eclipse/EclipseState/TracerConfig.hpp
diff --git a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
index a89ed3242..a2ee8aa84 100644
--- a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
+++ b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
@@ -22,6 +22,7 @@
 #define OPM_PARSER_ECLIPSE_GRID_HPP
 
 
+#include <opm/parser/eclipse/EclipseState/Grid/MapAxes.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/MinpvMode.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/PinchMode.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/GridDims.hpp>
@@ -206,9 +207,8 @@ namespace Opm {
         // Input grid data.
         std::vector<double> m_zcorn;
         std::vector<double> m_coord;
-        std::vector<double> m_mapaxes;
         std::vector<int> m_actnum;
-        std::string m_mapunits;
+        std::optional<MapAxes> m_mapaxes;
 
         // Mapping to/from active cells.
         int m_nactive;
diff --git a/opm/parser/eclipse/EclipseState/Grid/MapAxes.hpp b/opm/parser/eclipse/EclipseState/Grid/MapAxes.hpp
new file mode 100644
index 000000000..61ddcca32
--- /dev/null
+++ b/opm/parser/eclipse/EclipseState/Grid/MapAxes.hpp
@@ -0,0 +1,71 @@
+/*
+  Copyright 2021 Equinor ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef OPM_MAPAXES_HPP
+#define OPM_MAPAXES_HPP
+
+#include <array>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace Opm {
+
+/*
+  This class will internalize the information needed to transform (X,Y)
+  coordinates between grid input coordinates and world coordinates based on the
+  MAPAXES keyword.
+*/
+class Deck;
+
+namespace EclIO {
+class EclFile;
+}
+
+class MapAxes {
+public:
+    MapAxes();
+    MapAxes(double X1, double Y1, double X2, double Y2, double X3, double Y3);
+    MapAxes(const std::string& mapunits, double X1, double Y1, double X2, double Y2, double X3, double Y3);
+    MapAxes(EclIO::EclFile& egridfile);
+    MapAxes(const Deck& deck);
+
+    void transform(double& x, double& y) const;
+    void inv_transform(double& x, double& y) const;
+    const std::optional<std::string>& mapunits() const;
+    const std::vector<float>& input() const;
+    bool operator==(const MapAxes& other) const;
+
+private:
+    MapAxes(double length_factor, double X1, double Y1, double X2, double Y2, double X3, double Y3);
+    void init(double length_factor, double X1, double Y1, double X2, double Y2, double X3, double Y3);
+
+    std::array<double, 2> origin;
+    std::array<double, 2> unit_x;
+    std::array<double, 2> unit_y;
+    std::vector<float> m_input;
+    double inv_norm;
+    std::optional<std::string> map_units;
+};
+
+
+}
+
+#endif
diff --git a/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp b/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp
index 6c455bab5..87603976d 100644
--- a/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp
+++ b/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp
@@ -228,34 +228,8 @@ EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum)
 
     initGrid(deck);
 
-    if (deck.hasKeyword("MAPUNITS")){
-       if ((m_mapunits =="") || ( !keywInputBeforeGdfile(deck, "MAPUNITS"))) {
-            const auto& record = deck.getKeyword<ParserKeywords::MAPUNITS>( ).getStringData();
-            m_mapunits = record[0];
-       }
-    }
-
-    if (deck.hasKeyword("MAPAXES")){
-       if ((m_mapaxes.size() == 0) || ( !keywInputBeforeGdfile(deck, "MAPAXES"))) {
-          const Opm::DeckKeyword& mapaxesKeyword = deck.getKeyword<ParserKeywords::MAPAXES>();
-          const auto& mapaxes_data = mapaxesKeyword.getRawDoubleData();
-          if (mapaxes_data.size( ) != 6)
-              throw std::logic_error("Incorrect size for MAPAXES keyword");
-
-          auto length_factor = deck.getActiveUnitSystem().getDimension(UnitSystem::measure::length).getSIScaling();
-          if (m_mapunits == "METRES")
-              length_factor = unit::meter;
-          else if (m_mapunits == "FEET")
-              length_factor = unit::feet;
-          else if (m_mapunits == "CM")
-              length_factor = prefix::centi*unit::meter;
-
-          m_mapaxes.resize(6);
-          for (std::size_t n=0; n < 6; n++){
-              m_mapaxes[n] = length_factor * mapaxes_data[n];
-          }
-       }
-    }
+    if (deck.hasKeyword<ParserKeywords::MAPAXES>())
+        this->m_mapaxes = std::make_optional<MapAxes>( deck );
 
     if (actnum != nullptr) {
         resetACTNUM(actnum);
@@ -399,15 +373,9 @@ EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum)
             resetACTNUM( actnum );
         }
 
-        if (egridfile.hasKey("MAPAXES")) {
-            const std::vector<float>& mapaxes_f = egridfile.get<float>("MAPAXES");
-            m_mapaxes.assign(mapaxes_f.begin(), mapaxes_f.end());
-        }
+        if (egridfile.hasKey("MAPAXES"))
+            this->m_mapaxes = std::make_optional<MapAxes>(egridfile);
 
-        if (egridfile.hasKey("MAPUNITS")) {
-            const std::vector<std::string>& mapunits = egridfile.get<std::string>("MAPUNITS");
-            m_mapunits=mapunits[0];
-        }
 
         ZcornMapper mapper( getNX(), getNY(), getNZ());
         zcorn_fixed = mapper.fixupZCORN( m_zcorn );
@@ -1354,7 +1322,7 @@ std::vector<double> EclipseGrid::createDVector(const std::array<int,3>& dims, st
         if (m_zcorn.size() != other.m_zcorn.size())
             return false;
 
-        if (m_mapaxes.size() != other.m_mapaxes.size())
+        if (!(m_mapaxes == other.m_mapaxes))
             return false;
 
         if (m_actnum != other.m_actnum)
@@ -1366,9 +1334,6 @@ std::vector<double> EclipseGrid::createDVector(const std::array<int,3>& dims, st
         if (m_zcorn != other.m_zcorn)
             return false;
 
-        if (m_mapaxes != other.m_mapaxes)
-            return false;
-
         bool status = ((m_pinch == other.m_pinch)  && (m_minpvMode == other.getMinpvMode()));
 
         if(m_minpvMode!=MinpvMode::ModeEnum::Inactive) {
@@ -1641,8 +1606,6 @@ std::vector<double> EclipseGrid::createDVector(const std::array<int,3>& dims, st
             zcorn_f[n] = static_cast<double>(units.from_si(length, m_zcorn[n]));
         }
 
-        std::vector<float> mapaxes_f;
-
         std::vector<int> filehead(100,0);
         filehead[0] = 3;                     // version number
         filehead[1] = 2007;                  // release year
@@ -1682,32 +1645,8 @@ std::vector<double> EclipseGrid::createDVector(const std::array<int,3>& dims, st
             OPM_THROW(std::runtime_error, "Unit system not supported when writing to EGRID file");
             break;
         }
-
         gridunits.push_back("");
 
-        // map units is not dependent on deck units. A user may specify FIELD units for the model
-        // and metric units for the MAPAXES keyword (MAPUNITS)
-
-        std::vector<std::string> mapunits;
-
-        if ((m_mapunits.size() > 0) && (m_mapaxes.size() > 0)) {
-            mapunits.push_back(m_mapunits);
-        }
-
-        auto length_factor = units.getDimension(UnitSystem::measure::length).getSIScaling();
-        if (m_mapunits == "METRES")
-            length_factor = unit::meter;
-        else if (m_mapunits == "FEET")
-            length_factor = unit::feet;
-        else if (m_mapunits == "CM")
-            length_factor = prefix::centi*unit::meter;
-
-        if (m_mapaxes.size() > 0){
-            for (double dv :  m_mapaxes){
-                mapaxes_f.push_back(static_cast<float>(dv / length_factor));
-            }
-        }
-
         std::vector<int> endgrid = {};
 
         // Writing vectors to egrid file
@@ -1715,12 +1654,12 @@ std::vector<double> EclipseGrid::createDVector(const std::array<int,3>& dims, st
         Opm::EclIO::EclOutput egridfile(filename, formatted);
         egridfile.write("FILEHEAD", filehead);
 
-        if (mapunits.size() > 0) {
-            egridfile.write("MAPUNITS", mapunits);
-        }
+        if (this->m_mapaxes.has_value()) {
+            const auto& mapunits = this->m_mapaxes.value().mapunits();
+            if (mapunits.has_value())
+                egridfile.write("MAPUNITS", std::vector<std::string>{ mapunits.value() });
 
-        if (mapaxes_f.size() > 0){
-            egridfile.write("MAPAXES", mapaxes_f);
+            egridfile.write("MAPAXES", this->m_mapaxes.value().input());
         }
 
         egridfile.write("GRIDUNIT", gridunits);
diff --git a/src/opm/parser/eclipse/EclipseState/Grid/MapAxes.cpp b/src/opm/parser/eclipse/EclipseState/Grid/MapAxes.cpp
new file mode 100644
index 000000000..7fddabecf
--- /dev/null
+++ b/src/opm/parser/eclipse/EclipseState/Grid/MapAxes.cpp
@@ -0,0 +1,157 @@
+/*
+  Copyright 2021 Equinor ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <cmath>
+#include <stdexcept>
+#include <string>
+
+#include <opm/common/utility/String.hpp>
+#include <opm/parser/eclipse/Units/Units.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/MapAxes.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/io/eclipse/EclFile.hpp>
+#include <opm/io/eclipse/EclOutput.hpp>
+#include <opm/parser/eclipse/Parser/ParserKeywords/M.hpp>
+
+namespace Opm {
+
+namespace {
+
+double length_factor(const std::string& mapunits) {
+    auto mp = trim_copy(mapunits);
+    if (mp == "METRES")
+        return unit::meter;
+    else if (mp == "FEET")
+        return unit::feet;
+    else if (mp == "CM")
+        return prefix::centi*unit::meter;
+    throw std::logic_error("Unrecognized MAPUNITS value: " + mapunits);
+}
+
+}
+
+
+
+MapAxes::MapAxes() :
+    origin({0,0}),
+    unit_x({1,0}),
+    unit_y({0,1}),
+    m_input({1,0,0,0,0,1}),
+    inv_norm(1.0)
+{
+}
+
+
+const std::optional<std::string>& MapAxes::mapunits() const {
+    return this->map_units;
+}
+
+
+const std::vector<float>& MapAxes::input() const {
+    return this->m_input;
+}
+
+
+void MapAxes::init(double length_factor, double X1, double Y1, double X2, double Y2, double X3, double Y3)
+{
+    auto f = [](double d) { return static_cast<float>(d); };
+    this->m_input = {f(X1),f(Y1),f(X2),f(Y2),f(X3),f(Y3)};
+    this->origin = { length_factor*X2, length_factor*Y2};
+    this->unit_x = {X3 - X2, Y3 - Y2};
+    this->unit_y = {X1 - X2, Y1 - Y2};
+
+    auto norm_x = 1.0/std::hypot( this->unit_x[0], this->unit_x[1] );
+    auto norm_y = 1.0/std::hypot( this->unit_y[0], this->unit_y[1] );
+
+    this->unit_x[0] *= norm_x;
+    this->unit_x[1] *= norm_x;
+    this->unit_y[0] *= norm_y;
+    this->unit_y[1] *= norm_y;
+
+    this->inv_norm = 1.0 / (this->unit_x[0] * unit_y[1] - this->unit_x[1]*this->unit_y[0]);
+}
+
+MapAxes::MapAxes(double length_factor, double X1, double Y1, double X2, double Y2, double X3, double Y3) {
+    this->init(length_factor,X1,Y1,X2,Y2,X3,Y3);
+}
+
+MapAxes::MapAxes(double X1, double Y1, double X2, double Y2, double X3, double Y3):
+    MapAxes(1.0, X1, Y1, X2, Y2, X3, Y3)
+{}
+
+MapAxes::MapAxes(const std::string& mapunits, double X1, double Y1, double X2, double Y2, double X3, double Y3) :
+    MapAxes( length_factor(mapunits), X1, Y1, X2, Y2, X3, Y3)
+{
+    this->map_units = mapunits;
+}
+
+
+MapAxes::MapAxes(const Deck& deck)
+{
+    if (!deck.hasKeyword<ParserKeywords::MAPAXES>())
+        throw std::logic_error("Can not instantiate MapAxes object without MAPAXES keyword in deck");
+
+    const auto& mapaxes_kw = deck.getKeyword<ParserKeywords::MAPAXES>();
+    const auto& mapaxes_data = mapaxes_kw.getRawDoubleData();
+    double lf = 1.0;
+    if (deck.hasKeyword<ParserKeywords::MAPUNITS>()) {
+        this->map_units = deck.getKeyword<ParserKeywords::MAPUNITS>().getRecord(0).getItem(0).get<std::string>(0);
+        lf = length_factor(this->map_units.value());
+    }
+
+    this->init(lf, mapaxes_data[0], mapaxes_data[1], mapaxes_data[2], mapaxes_data[3], mapaxes_data[4], mapaxes_data[5]);
+}
+
+MapAxes::MapAxes(EclIO::EclFile& egridfile)
+{
+    if (!egridfile.hasKey("MAPAXES"))
+        throw std::logic_error("Can not instantiate MapAxes object without MAPAXES keyword in EGRID");
+
+    const auto& mapaxes_data = egridfile.get<float>("MAPAXES");
+    double lf = 1.0;
+    if (egridfile.hasKey("MAPUNITS")) {
+        this->map_units = egridfile.get<std::string>("MAPUNITS")[0];
+        lf = length_factor(this->map_units.value());
+    }
+    this->init(lf, mapaxes_data[0], mapaxes_data[1], mapaxes_data[2], mapaxes_data[3], mapaxes_data[4], mapaxes_data[5]);
+ }
+
+
+void MapAxes::transform(double& x, double& y) const {
+    double tmpx = x;
+    x = this->origin[0] + tmpx*this->unit_x[0] + y*this->unit_y[0];
+    y = this->origin[1] + tmpx*this->unit_x[1] + y*this->unit_y[1];
+}
+
+
+void MapAxes::inv_transform(double& x, double& y) const {
+    auto dx = x - this->origin[0];
+    auto dy = y - this->origin[1];
+
+    x = this->inv_norm * ( dx * this->unit_y[1] - dy*this->unit_y[0]);
+    y = this->inv_norm * (-dx * this->unit_x[1] + dy*this->unit_x[0]);
+}
+
+bool MapAxes::operator==(const MapAxes& other) const {
+    return this->m_input == other.m_input &&
+           this->map_units == other.map_units;
+}
+
+}
+
diff --git a/tests/parser/EclipseGridTests.cpp b/tests/parser/EclipseGridTests.cpp
index ebc262c1d..f9540a9aa 100644
--- a/tests/parser/EclipseGridTests.cpp
+++ b/tests/parser/EclipseGridTests.cpp
@@ -29,6 +29,7 @@
 
 #include <opm/common/utility/FileSystem.hpp>
 #include <opm/common/utility/OpmInputError.hpp>
+#include <opm/parser/eclipse/Units/Units.hpp>
 
 #define BOOST_TEST_MODULE EclipseGridTests
 #include <boost/test/unit_test.hpp>
@@ -41,6 +42,7 @@
 
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/MapAxes.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/GridDims.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/PinchMode.hpp>
@@ -2877,3 +2879,72 @@ BOOST_AUTO_TEST_CASE(TEST_COLLAPSED_CELL) {
     for (std::size_t g = 0; g < grid.getCartesianSize(); g++)
         BOOST_CHECK_EQUAL(grid.getCellVolume(g), 0);
 }
+
+
+
+BOOST_AUTO_TEST_CASE(MAPAXES_DEFAULT) {
+    Opm::MapAxes ma1;
+    const double x1 = 77;
+    const double y1 = 99;
+    {
+        double xt{x1};
+        double yt{y1};
+
+        ma1.transform(xt,yt);
+        BOOST_CHECK_EQUAL(x1,xt);
+        BOOST_CHECK_EQUAL(y1,yt);
+    }
+    {
+        double xi{x1};
+        double yi{y1};
+
+        ma1.inv_transform(xi,yi);
+        BOOST_CHECK_EQUAL(x1,xi);
+        BOOST_CHECK_EQUAL(y1,yi);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(MAPAXES) {
+    const double origo_x = 10;
+    const double origo_y = 10;
+    const double px = 5;
+    const double py = 5;
+    Opm::MapAxes ma( 5,5, origo_x, origo_y, 15,5 );
+    double xt{px};
+    double yt{py};
+    ma.transform(xt,yt);
+    BOOST_CHECK_EQUAL( xt, origo_x );
+    BOOST_CHECK_CLOSE( yt, origo_y - std::sqrt(2) * px, 1e-8);
+
+    double xi{xt};
+    double yi{yt};
+    ma.inv_transform(xi,yi);
+    BOOST_CHECK_CLOSE( px, xi, 1e-8);
+    BOOST_CHECK_CLOSE( py, yi, 1e-8);
+}
+
+BOOST_AUTO_TEST_CASE(MAPAXES_MAPUNITS) {
+    const double feet = Opm::unit::feet;
+    const double origo_x = 10;
+    const double origo_y = 10;
+    const double px = 5;
+    const double py = 5;
+    Opm::MapAxes ma( "FEET", 5/feet,5/feet, origo_x/feet, origo_y/feet, 15/feet, 5/feet );
+    double xt{px};
+    double yt{py};
+    ma.transform(xt,yt);
+    BOOST_CHECK_CLOSE( xt, origo_x , 1e-8);
+    BOOST_CHECK_CLOSE( yt, origo_y - std::sqrt(2) * px, 1e-8);
+
+    double xi{xt};
+    double yi{yt};
+    ma.inv_transform(xi,yi);
+    BOOST_CHECK_CLOSE( px, xi, 1e-8);
+    BOOST_CHECK_CLOSE( py, yi, 1e-8);
+
+    auto mapunits = ma.mapunits();
+    BOOST_CHECK_EQUAL(mapunits.value(), "FEET");
+    BOOST_CHECK_THROW(Opm::MapAxes( "NO_SUCH_LENGTH_UNIT", 5,5, origo_x, origo_y, 15, 5 ), std::exception);
+
+}
+