diff --git a/src/mdarray/group.jl b/src/mdarray/group.jl
index 45d3c6b3..6d56e64e 100644
--- a/src/mdarray/group.jl
+++ b/src/mdarray/group.jl
@@ -235,7 +235,7 @@ function unsafe_createmdarray(
         group,
         name,
         length(dimensions),
-        DimensionHList(dimensions),
+        DimensionHList(reverse(dimensions)),
         datatype,
         CSLConstListWrapper(options),
     )
@@ -257,7 +257,7 @@ function createmdarray(
         group,
         name,
         length(dimensions),
-        DimensionHList(dimensions),
+        DimensionHList(reverse(dimensions)),
         datatype,
         CSLConstListWrapper(options),
     )
diff --git a/src/mdarray/highlevel.jl b/src/mdarray/highlevel.jl
index eb7be715..55ae1869 100644
--- a/src/mdarray/highlevel.jl
+++ b/src/mdarray/highlevel.jl
@@ -7,8 +7,7 @@ function writemdarray(
     options::OptionList = nothing,
 )::Nothing where {T<:NumericAttributeType,D}
     dimensions = AbstractDimension[
-        createdimension(group, "$name.$d", "", "", size(value, d)) for
-        d in D:-1:1
+        createdimension(group, "$name.$d", "", "", size(value, d)) for d in 1:D
     ]
     extendeddatatypecreate(T) do datatype
         createmdarray(group, name, dimensions, datatype, options) do mdarray
diff --git a/src/mdarray/mdarray.jl b/src/mdarray/mdarray.jl
index f5ea4556..a5ccb9d1 100644
--- a/src/mdarray/mdarray.jl
+++ b/src/mdarray/mdarray.jl
@@ -254,7 +254,7 @@ function unsafe_getindex(
     indices::Integer...,
 )::AbstractMDArray
     @assert !isnull(mdarray)
-    viewexpr = "[" * join(indices, ",") * "]"
+    viewexpr = "[" * join(reverse(indices), ",") * "]"
     return unsafe_getview(mdarray, viewexpr)
 end
 
@@ -263,7 +263,7 @@ function getindex(
     indices::Integer...,
 )::AbstractMDArray
     @assert !isnull(mdarray)
-    viewexpr = "[" * join(indices, ",") * "]"
+    viewexpr = "[" * join(reverse(indices), ",") * "]"
     return getview(mdarray, viewexpr)
 end
 
@@ -348,7 +348,7 @@ function unsafe_getresampled(
     ptr = GDAL.gdalmdarraygetresampled(
         mdarray,
         isnothing(newdims) ? 0 : length(newdims),
-        isnothing(newdims) ? C_NULL : DimensionHList(newdims),
+        isnothing(newdims) ? C_NULL : DimensionHList(reverse(newdims)),
         resamplealg,
         isnothing(targetsrs) ? C_NULL : targetsrs,
         CSLConstListWrapper(options),
@@ -368,7 +368,7 @@ function getresampled(
     ptr = GDAL.gdalmdarraygetresampled(
         mdarray,
         isnothing(newdims) ? 0 : length(newdims),
-        isnothing(newdims) ? C_NULL : DimensionHList(newdims),
+        isnothing(newdims) ? C_NULL : DimensionHList(reverse(newdims)),
         resamplealg,
         isnothing(targetsrs) ? C_NULL : targetsrs,
         CSLConstListWrapper(options),
@@ -563,7 +563,7 @@ function adviseread(
     @assert !isnull(mdarray)
     return GDAL.gdalmdarrayadviseread(
         mdarray,
-        isnothing(arraystartidx) ? C_NULL : arraystartidx,
+        isnothing(arraystartidx) ? C_NULL : reverse(arraystartidx),
         isnothing(count) ? C_NULL : count,
         CSLConstListWrapper(options),
     )
@@ -647,7 +647,7 @@ function getdimensions(
     dimensionshptr = GDAL.gdalmdarraygetdimensions(mdarray, dimensionscountref)
     dimensions = AbstractDimension[
         IDimension(unsafe_load(dimensionshptr, n), mdarray.dataset) for
-        n in 1:dimensionscountref[]
+        n in dimensionscountref[]:-1:1
     ]
     GDAL.vsifree(dimensionshptr)
     return dimensions
@@ -661,7 +661,7 @@ function unsafe_getdimensions(
     dimensionshptr = GDAL.gdalmdarraygetdimensions(mdarray, dimensionscountref)
     dimensions = AbstractDimension[
         Dimension(unsafe_load(dimensionshptr, n), mdarray.dataset) for
-        n in 1:dimensionscountref[]
+        n in dimensionscountref[]:-1:1
     ]
     GDAL.vsifree(dimensionshptr)
     return dimensions
@@ -670,7 +670,7 @@ end
 function Base.size(mdarray::AbstractMDArray)
     getdimensions(mdarray) do dimensions
         D = length(dimensions)
-        return ntuple(d -> getsize(dimensions[D-d+1]), D)
+        return ntuple(d -> getsize(dimensions[d]), D)
     end
 end
 
@@ -800,7 +800,7 @@ end
 function read(mdarray::AbstractMDArray)::AbstractArray
     getdimensions(mdarray) do dimensions
         D = length(dimensions)
-        sz = [getsize(dimensions[d]) for d in D:-1:1]
+        sz = [getsize(dimensions[d]) for d in 1:D]
         getdatatype(mdarray) do datatype
             class = getclass(datatype)
             @assert class == GDAL.GEDTC_NUMERIC
diff --git a/test/test_mdarray.jl b/test/test_mdarray.jl
index 12347d95..102240d5 100644
--- a/test/test_mdarray.jl
+++ b/test/test_mdarray.jl
@@ -215,7 +215,7 @@ end
 
                 # @test AG.iswritable(mdarray)
 
-                data = Float32[100 * x + y for y in 1:ny, x in 1:nx]
+                data = Float32[x + 100 * y for x in 1:nx, y in 1:ny]
                 AG.write(mdarray, data)
 
                 write_attributes(mdarray)
@@ -308,12 +308,12 @@ end
                 @test AG.getclass(datatype) == GDAL.GEDTC_NUMERIC
                 @test AG.getnumericdatatype(datatype) == AG.GDT_Float32
 
-                data = Array{Float32}(undef, ny, nx)
+                data = Array{Float32}(undef, nx, ny)
                 AG.read!(mdarray, data)
-                @test data == Float32[100 * x + y for y in 1:ny, x in 1:nx]
+                @test data == Float32[x + 100 * y for x in 1:nx, y in 1:ny]
 
                 data = AG.read(mdarray)
-                @test data == Float32[100 * x + y for y in 1:ny, x in 1:nx]
+                @test data == Float32[x + 100 * y for x in 1:nx, y in 1:ny]
 
                 test_attributes(mdarray)
 
@@ -411,8 +411,8 @@ end
                                             # @test AG.iswritable(mdarray)
 
                                             data = Float32[
-                                                100 * x + y for y in 1:ny,
-                                                x in 1:nx
+                                                x + 100 * y for x in 1:nx,
+                                                y in 1:ny
                                             ]
                                             AG.write(mdarray, data)
 
@@ -539,15 +539,15 @@ end
                                         ) do doesnotexist
                                         end
 
-                                        data = Array{Float32}(undef, ny, nx)
+                                        data = Array{Float32}(undef, nx, ny)
                                         AG.read!(mdarray, data)
                                         @test data == Float32[
-                                            100 * x + y for y in 1:ny, x in 1:nx
+                                            x + 100 * y for x in 1:nx, y in 1:ny
                                         ]
 
                                         data = AG.read(mdarray)
                                         @test data == Float32[
-                                            100 * x + y for y in 1:ny, x in 1:nx
+                                            x + 100 * y for x in 1:nx, y in 1:ny
                                         ]
 
                                         test_attributes(mdarray)