diff --git a/autotest/ogr/ogr_pg.py b/autotest/ogr/ogr_pg.py
index d1cc694faae9..a43f989870d4 100755
--- a/autotest/ogr/ogr_pg.py
+++ b/autotest/ogr/ogr_pg.py
@@ -1526,7 +1526,7 @@ def test_ogr_pg_32():
     # Create third layer with very approximative EPSG:4326 but without authority
 
     srs = osr.SpatialReference()
-    srs.SetFromUserInput("""GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]""")
+    srs.SetFromUserInput("""GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]""")
     gdaltest.pg_lyr = gdaltest.pg_ds.CreateLayer('testsrtext3', srs=srs)
 
     # Must still be 1
diff --git a/autotest/ogr/ogr_sqlite.py b/autotest/ogr/ogr_sqlite.py
index d13f7f1aa6b5..ca6a0611e64c 100755
--- a/autotest/ogr/ogr_sqlite.py
+++ b/autotest/ogr/ogr_sqlite.py
@@ -818,7 +818,7 @@ def test_ogr_sqlite_17(require_spatialite):
     assert lyr is None, 'layer creation should have failed'
 
     srs = osr.SpatialReference()
-    srs.SetFromUserInput("""GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]""")
+    srs.ImportFromEPSG(4326)
     lyr = ds.CreateLayer('geomspatialite', srs=srs)
 
     geom = ogr.CreateGeometryFromWkt('POINT(0 1)')
diff --git a/autotest/osr/osr_epsg.py b/autotest/osr/osr_epsg.py
index 60314c6d4c70..d74b58ff0dec 100755
--- a/autotest/osr/osr_epsg.py
+++ b/autotest/osr/osr_epsg.py
@@ -32,7 +32,7 @@
 
 
 
-from osgeo import osr
+from osgeo import ogr, osr
 import gdaltest
 import pytest
 
@@ -445,3 +445,14 @@ def test_osr_epsg_auto_identify_epsg_odd_wkt():
     srs = osr.SpatialReference()
     srs.SetFromUserInput("""PROJCS["GDA94 / MGA zone 56 / AUSGeoid09_GDA94_V1",GEOGCS["GDA94 / MGA zone 56 / AUSGeoid09_GDA94_V1",DATUM["GDA94",SPHEROID["GRS 1980",6378137.000,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0.0000000000000000,AUTHORITY["EPSG","8901"]],UNIT["Degree",0.01745329251994329547,AUTHORITY["EPSG","9102"]],AUTHORITY["EPSG","28356"]],UNIT["Meter",1.00000000000000000000,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0.0000000000000000],PARAMETER["central_meridian",153.0000000000000000],PARAMETER["scale_factor",0.9996000000000000],PARAMETER["false_easting",500000.000],PARAMETER["false_northing",10000000.000],AXIS["Easting",EAST],AXIS["Northing",NORTH],AXIS["Height",UP],AUTHORITY["EPSG","28356"]]""")
     srs.AutoIdentifyEPSG()
+
+###############################################################################
+#   Test bugfix for https://github.com/OSGeo/gdal/issues/4038
+
+
+def test_osr_epsg_auto_identify_epsg_projcrs_with_geogcrs_without_axis_roder():
+
+    srs = osr.SpatialReference()
+    assert srs.SetFromUserInput("""PROJCS["GDA_1994_MGA_Zone_55",GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",147.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]""") == 0
+    assert srs.AutoIdentifyEPSG() != ogr.OGRERR_NONE
+    assert srs.CloneGeogCS().GetAuthorityCode(None) is None
diff --git a/gdal/ogr/ogrspatialreference.cpp b/gdal/ogr/ogrspatialreference.cpp
index 7d389bb2d273..afd71e147b14 100644
--- a/gdal/ogr/ogrspatialreference.cpp
+++ b/gdal/ogr/ogrspatialreference.cpp
@@ -11692,11 +11692,40 @@ OGRErr OSRDemoteTo2D( OGRSpatialReferenceH hSRS, const char* pszName  )
 int OGRSpatialReference::GetEPSGGeogCS() const
 
 {
-    const char *pszAuthName = GetAuthorityName( "GEOGCS" );
+/* -------------------------------------------------------------------- */
+/*      Check axis order.                                               */
+/* -------------------------------------------------------------------- */
+    auto poGeogCRS = std::unique_ptr<OGRSpatialReference>(CloneGeogCS());
+    if( !poGeogCRS )
+        return -1;
+
+    bool ret = false;
+    poGeogCRS->d->demoteFromBoundCRS();
+    auto cs = proj_crs_get_coordinate_system(d->getPROJContext(),
+                                             poGeogCRS->d->m_pj_crs);
+    poGeogCRS->d->undoDemoteFromBoundCRS();
+    if( cs )
+    {
+        const char* pszDirection = nullptr;
+        if( proj_cs_get_axis_info(
+            d->getPROJContext(), cs, 0, nullptr, nullptr, &pszDirection,
+            nullptr, nullptr, nullptr, nullptr) )
+        {
+            if( EQUAL(pszDirection, "north") )
+            {
+                ret = true;
+            }
+        }
+
+        proj_destroy(cs);
+    }
+    if( !ret )
+        return -1;
 
 /* -------------------------------------------------------------------- */
 /*      Do we already have it?                                          */
 /* -------------------------------------------------------------------- */
+    const char *pszAuthName = GetAuthorityName( "GEOGCS" );
     if( pszAuthName != nullptr && EQUAL(pszAuthName, "epsg") )
         return atoi(GetAuthorityCode( "GEOGCS" ));