diff --git a/src/gdal.cpp b/src/gdal.cpp
index 5f6cf9b90..ad92199de 100644
--- a/src/gdal.cpp
+++ b/src/gdal.cpp
@@ -17,6 +17,16 @@
 #include "wkb.h"
 #include "gdal_sf_pkg.h"
 
+#if GDAL_VERSION_MAJOR == 2
+# if GDAL_VERSION_MINOR >= 5
+#  define HAVE250
+# endif
+#else
+# if GDAL_VERSION_MAJOR > 2
+#  define HAVE250
+# endif
+#endif
+
 //
 // Returns errors to R
 // Note only case 4 actually returns immediately
@@ -111,6 +121,7 @@ Rcpp::List CPL_crs_parameters(std::string p4s) {
 	Rcpp::List out(7);
 	OGRErr Err;
 	OGRSpatialReference *srs = new OGRSpatialReference;
+	srs = handle_axis_order(srs);
 	handle_error(srs->importFromProj4(p4s.c_str()));
 	out(0) = Rcpp::NumericVector::create(srs->GetSemiMajor());
 	out(1) = Rcpp::NumericVector::create(srs->GetSemiMinor());
@@ -138,8 +149,10 @@ Rcpp::List CPL_crs_parameters(std::string p4s) {
 Rcpp::LogicalVector CPL_crs_equivalent(std::string crs1, std::string crs2) {
 	Rcpp::LogicalVector out(1);
 	OGRSpatialReference *srs1 = new OGRSpatialReference;
+	srs1 = handle_axis_order(srs1);
 	handle_error(srs1->importFromProj4(crs1.c_str()));
 	OGRSpatialReference *srs2 = new OGRSpatialReference;
+	srs2 = handle_axis_order(srs2);
 	handle_error(srs2->importFromProj4(crs2.c_str()));
 	out(0) = (bool) srs1->IsSame(srs2);
 	delete srs1;
@@ -158,6 +171,7 @@ std::vector<OGRGeometry *> ogr_from_sfc(Rcpp::List sfc, OGRSpatialReference **sr
 	if (p4s != NA_STRING) {
 		Rcpp::CharacterVector cv = crs["proj4string"];
 		local_srs = new OGRSpatialReference;
+		local_srs = handle_axis_order(local_srs);
 		OGRErr err = local_srs->importFromProj4(cv[0]);
 		if (err != OGRERR_NONE) {
 			local_srs->Release(); // #nocov
@@ -272,6 +286,7 @@ Rcpp::List sfc_from_ogr(std::vector<OGRGeometry *> g, bool destroy = false) {
 // [[Rcpp::export]]
 Rcpp::List CPL_crs_from_epsg(int epsg) {
 	OGRSpatialReference ref;
+	handle_axis_order(&ref);
 	if (ref.importFromEPSG(epsg) == OGRERR_NONE)
 		return get_crs(&ref);
 	else
@@ -282,6 +297,7 @@ Rcpp::List CPL_crs_from_epsg(int epsg) {
 Rcpp::List CPL_crs_from_wkt(Rcpp::CharacterVector wkt) {
 	char *cp = wkt[0];
 	OGRSpatialReference ref;
+	handle_axis_order(&ref);
 #if GDAL_VERSION_MAJOR <= 2 && GDAL_VERSION_MINOR <= 2
 	handle_error(ref.importFromWkt(&cp));
 #else
@@ -359,6 +375,7 @@ Rcpp::List CPL_transform(Rcpp::List sfc, Rcpp::CharacterVector proj4) {
 
 	// import proj4string:
 	OGRSpatialReference *dest = new OGRSpatialReference;
+	dest = handle_axis_order(dest);
 	handle_error(dest->importFromProj4((const char *) (proj4[0])));
 
 	// transform geometries:
@@ -410,6 +427,7 @@ Rcpp::List CPL_wrap_dateline(Rcpp::List sfc, Rcpp::CharacterVector opt, bool qui
 // [[Rcpp::export]]
 Rcpp::List CPL_crs_from_proj4string(Rcpp::CharacterVector p4s) {
 	OGRSpatialReference ref;
+	handle_axis_order(&ref);
 	if (ref.importFromProj4(p4s[0]) == OGRERR_NONE)
 		return get_crs(&ref);
 	else {
@@ -470,3 +488,15 @@ Rcpp::LogicalVector CPL_gdal_with_geos() {
 	bool withGEOS = OGRGeometryFactory::haveGEOS(); 
 	return Rcpp::LogicalVector::create(withGEOS);
 }
+
+OGRSpatialReference *handle_axis_order(OGRSpatialReference *sr, bool traditional) {
+#ifdef HAVE250
+	if (sr != NULL) {
+		if (traditional)
+			sr->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+		else
+			sr->SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT);
+	}
+#endif
+	return sr;
+}
diff --git a/src/gdal.h b/src/gdal.h
index 7afedacc3..d7df9412f 100644
--- a/src/gdal.h
+++ b/src/gdal.h
@@ -1,2 +1,3 @@
 void set_error_handler(void);
 void unset_error_handler(void);
+OGRSpatialReference *handle_axis_order(OGRSpatialReference *sr, bool traditional = true);
diff --git a/src/gdal_write.cpp b/src/gdal_write.cpp
index 722821cc4..55150bcb6 100644
--- a/src/gdal_write.cpp
+++ b/src/gdal_write.cpp
@@ -246,6 +246,7 @@ int CPL_write_ogr(Rcpp::List obj, Rcpp::CharacterVector dsn, Rcpp::CharacterVect
 	// read geometries:
 	OGRSpatialReference *sref = NULL;
 	std::vector<OGRGeometry *> geomv = ogr_from_sfc(geom, &sref);
+	sref = handle_axis_order(sref);
 
 	// create layer:
 	options = create_options(lco, quiet);
diff --git a/src/polygonize.cpp b/src/polygonize.cpp
index 917068965..16b79097a 100644
--- a/src/polygonize.cpp
+++ b/src/polygonize.cpp
@@ -76,6 +76,7 @@ Rcpp::List CPL_polygonize(Rcpp::CharacterVector raster, Rcpp::CharacterVector ma
 		Rcpp::stop("Creation failed.\n"); // #nocov
 	}
 	OGRSpatialReference *sr = new OGRSpatialReference;
+	sr = handle_axis_order(sr);
 	char **ppt = (char **) &wkt;
 #if GDAL_VERSION_MAJOR <= 2 && GDAL_VERSION_MINOR <= 2
 	sr->importFromWkt(ppt);
diff --git a/src/stars.cpp b/src/stars.cpp
index 37851bd4e..35e9e9bc6 100644
--- a/src/stars.cpp
+++ b/src/stars.cpp
@@ -8,6 +8,7 @@
 
 #include <Rcpp.h>
 
+#include "gdal.h"
 #include "gdal_read.h"
 #include "gdal_sf_pkg.h"
 
@@ -229,6 +230,7 @@ List CPL_read_gdal(CharacterVector fname, CharacterVector options, CharacterVect
 		proj = CharacterVector::create(wkt);
 		// proj4string:
 		OGRSpatialReference *sr = new OGRSpatialReference;
+		sr = handle_axis_order(sr);
 		char **ppt = (char **) &wkt;
 #if GDAL_VERSION_MAJOR <= 2 && GDAL_VERSION_MINOR <= 2
 		sr->importFromWkt(ppt);