From 64c0a48b80a6d8659468ca5771d46ecceb76b916 Mon Sep 17 00:00:00 2001 From: Fabian Gans Date: Thu, 21 May 2026 10:43:58 +0200 Subject: [PATCH 1/6] Only revert xy for epsg 4326 --- ext/ArchGDALExt/ArchGDALExt.jl | 95 ++++++++++++++++++++-------------- 1 file changed, 55 insertions(+), 40 deletions(-) diff --git a/ext/ArchGDALExt/ArchGDALExt.jl b/ext/ArchGDALExt/ArchGDALExt.jl index e21416d..7efeb91 100644 --- a/ext/ArchGDALExt/ArchGDALExt.jl +++ b/ext/ArchGDALExt/ArchGDALExt.jl @@ -2,7 +2,7 @@ module ArchGDALExt import ArchGDAL: RasterDataset, AbstractRasterBand, getgeotransform, width, height import ArchGDAL: getname, getcolorinterp, getband, nraster, getdataset using ArchGDAL: ArchGDAL as AG - +using ArchGDAL.GeoFormatTypes: WellKnownText, EPSG import ArchGDAL.DiskArrays: GridChunks, eachchunk import ArchGDAL.DiskArrays @@ -13,16 +13,20 @@ module ArchGDALExt include("archgdaldataset.jl") function dimname(a::RasterDataset, i) + xy = dimsymbols(a) if i == 1 - return :Y + return xy[1] elseif i == 2 - return :X + return xy[2] elseif i == 3 return :Band else error("RasterDataset only has 3 dimensions") end end + + + function dimvals(a::RasterDataset, i) if i == 1 geo=getgeotransform(a) @@ -41,6 +45,16 @@ module ArchGDALExt end end +function dimsymbols(a) + wkt = WellKnownText(AG.toWKT(AG.newspatialref(AG.getproj(a)))) + epsg = convert(EPSG, wkt) + if epsg.val == 4326 + (:Y, :X) + else + (:X, :Y) + end +end + iscontdim(a::RasterDataset, i) = i < 3 ? true : nraster(a)<8 function getattributes(a::RasterDataset) globatts = Dict{String,Any}( @@ -55,45 +69,46 @@ module ArchGDALExt end - function dimname(::AbstractRasterBand, i) - if i == 1 - return :Y - elseif i == 2 - return :X - else - error("RasterDataset only has 3 dimensions") - end +function dimname(a::AbstractRasterBand, i) + xy = dimsymbols(a) + if i == 1 + return xy[1] + elseif i == 2 + return xy[2] + else + error("RasterDataset only has 3 dimensions") end - function dimvals(b::AbstractRasterBand, i) - geo = getgeotransform(getdataset(b)) - if i == 1 - range(geo[1],length=width(b), step=geo[2]) - elseif i == 2 - range(geo[4],length=height(b), step=geo[6]) - else - error("RasterDataset only has 3 dimensions") - end - end - iscontdim(a::AbstractRasterBand, i) = true - function getattributes(a::AbstractRasterBand) - atts = getattributes(AG.RasterDataset(AG.getdataset(a))) - bandatts = getbandattributes(a) - return merge(atts, bandatts) +end +function dimvals(b::AbstractRasterBand, i) + geo = getgeotransform(getdataset(b)) + if i == 1 + range(geo[1], length=width(b), step=geo[2]) + elseif i == 2 + range(geo[4], length=height(b), step=geo[6]) + else + error("RasterDataset only has 3 dimensions") end +end +iscontdim(a::AbstractRasterBand, i) = true +function getattributes(a::AbstractRasterBand) + atts = getattributes(AG.RasterDataset(AG.getdataset(a))) + bandatts = getbandattributes(a) + return merge(atts, bandatts) +end - function insertattifnot!(attrs, val, name, condition) - if !condition(val) - attrs[name] = val - end - end - function getbandattributes(a::AbstractRasterBand) - atts = Dict{String,Any}() - catdict = Dict((i-1)=>v for (i,v) in enumerate(AG.getcategorynames(a))) - insertattifnot!(atts, AG.getnodatavalue(a), "missing_value", isnothing) - insertattifnot!(atts, catdict, "labels", isempty) - insertattifnot!(atts, AG.getunittype(a), "units", isempty) - insertattifnot!(atts, AG.getoffset(a), "add_offset", iszero) - insertattifnot!(atts, AG.getscale(a), "scale_factor", x->isequal(x, one(x))) - return atts +function insertattifnot!(attrs, val, name, condition) + if !condition(val) + attrs[name] = val end +end +function getbandattributes(a::AbstractRasterBand) + atts = Dict{String,Any}() + catdict = Dict((i - 1) => v for (i, v) in enumerate(AG.getcategorynames(a))) + insertattifnot!(atts, AG.getnodatavalue(a), "missing_value", isnothing) + insertattifnot!(atts, catdict, "labels", isempty) + insertattifnot!(atts, AG.getunittype(a), "units", isempty) + insertattifnot!(atts, AG.getoffset(a), "add_offset", iszero) + insertattifnot!(atts, AG.getscale(a), "scale_factor", x -> isequal(x, one(x))) + return atts +end end \ No newline at end of file From feb6fd52af190c538d7d08c3a12dfc956bf194bf Mon Sep 17 00:00:00 2001 From: Fabian Gans Date: Thu, 21 May 2026 15:17:11 +0200 Subject: [PATCH 2/6] call getproj only on datasets --- ext/ArchGDALExt/ArchGDALExt.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ext/ArchGDALExt/ArchGDALExt.jl b/ext/ArchGDALExt/ArchGDALExt.jl index 7efeb91..a8646d9 100644 --- a/ext/ArchGDALExt/ArchGDALExt.jl +++ b/ext/ArchGDALExt/ArchGDALExt.jl @@ -45,7 +45,8 @@ using ArchGDAL.GeoFormatTypes: WellKnownText, EPSG end end -function dimsymbols(a) +dimsymbols(a::AbstractRasterBand) = dimsymbols(AG.getdataset(a)) +function dimsymbols(a::RasterDataset) wkt = WellKnownText(AG.toWKT(AG.newspatialref(AG.getproj(a)))) epsg = convert(EPSG, wkt) if epsg.val == 4326 From b77ed0c743a95c219b0ca5279d0ef33b4cd084e1 Mon Sep 17 00:00:00 2001 From: Fabian Gans Date: Thu, 21 May 2026 15:21:26 +0200 Subject: [PATCH 3/6] fix some tests --- test/arrays.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/arrays.jl b/test/arrays.jl index 9794713..0cd12d5 100644 --- a/test/arrays.jl +++ b/test/arrays.jl @@ -62,9 +62,9 @@ end using ArchGDAL AG=ArchGDAL r = AG.readraster(p) - @test YAXArrayBase.dimnames(r) == (:Y, :X, :Band) - @test YAXArrayBase.dimname(r,1) == :Y - @test YAXArrayBase.dimname(r,2) == :X + @test YAXArrayBase.dimnames(r) == (:X, :Y, :Band) + @test YAXArrayBase.dimname(r, 1) == :X + @test YAXArrayBase.dimname(r, 2) == :Y @test YAXArrayBase.dimname(r,3) == :Band @test YAXArrayBase.dimvals(r,1) == -28493.166784412522:60.02213698319374:2298.189487965865 @test YAXArrayBase.dimvals(r,2) == 4.2558845438021915e6:-60.02213698319374:4.22503316539283e6 @@ -76,9 +76,9 @@ end @test YAXArrayBase.iscontdim(r,3) == true @test YAXArrayBase.getattributes(r)["projection_PROJ4"] == "+proj=cea +lat_ts=33.75 +lon_0=-117.333333333333 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs" b = AG.getband(r,1) - @test YAXArrayBase.dimnames(b) == (:Y, :X) - @test YAXArrayBase.dimname(b,1) == :Y - @test YAXArrayBase.dimname(b,2) == :X + @test YAXArrayBase.dimnames(b) == (:X, :Y) + @test YAXArrayBase.dimname(b, 1) == :X + @test YAXArrayBase.dimname(b, 2) == :Y @test YAXArrayBase.dimvals(b,1) == -28493.166784412522:60.02213698319374:2298.189487965865 @test YAXArrayBase.dimvals(b,2) == 4.2558845438021915e6:-60.02213698319374:4.22503316539283e6 @test_throws Exception YAXArrayBase.dimname(b,3) From 9b1cabc9d6084742b8a0f7f2fe7cada5dc8c7b12 Mon Sep 17 00:00:00 2001 From: Fabian Gans Date: Thu, 21 May 2026 15:27:22 +0200 Subject: [PATCH 4/6] one more try for fix --- ext/ArchGDALExt/ArchGDALExt.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ext/ArchGDALExt/ArchGDALExt.jl b/ext/ArchGDALExt/ArchGDALExt.jl index a8646d9..a40ad2e 100644 --- a/ext/ArchGDALExt/ArchGDALExt.jl +++ b/ext/ArchGDALExt/ArchGDALExt.jl @@ -46,8 +46,10 @@ using ArchGDAL.GeoFormatTypes: WellKnownText, EPSG end dimsymbols(a::AbstractRasterBand) = dimsymbols(AG.getdataset(a)) -function dimsymbols(a::RasterDataset) - wkt = WellKnownText(AG.toWKT(AG.newspatialref(AG.getproj(a)))) +get_projection(a) = AG.getproj(a) +get_projection(a::AbstractRasterBand) = AG.getproj(AG.getdataset(a)) +function dimsymbols(a) + wkt = WellKnownText(AG.toWKT(AG.newspatialref(get_projection(a)))) epsg = convert(EPSG, wkt) if epsg.val == 4326 (:Y, :X) From bc6fb4b6cc898fb4fde46838a6ba00c04e544fea Mon Sep 17 00:00:00 2001 From: Fabian Gans Date: Thu, 21 May 2026 16:36:58 +0200 Subject: [PATCH 5/6] remove special case again --- ext/ArchGDALExt/ArchGDALExt.jl | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/ext/ArchGDALExt/ArchGDALExt.jl b/ext/ArchGDALExt/ArchGDALExt.jl index a40ad2e..8f770d7 100644 --- a/ext/ArchGDALExt/ArchGDALExt.jl +++ b/ext/ArchGDALExt/ArchGDALExt.jl @@ -45,17 +45,21 @@ using ArchGDAL.GeoFormatTypes: WellKnownText, EPSG end end -dimsymbols(a::AbstractRasterBand) = dimsymbols(AG.getdataset(a)) -get_projection(a) = AG.getproj(a) -get_projection(a::AbstractRasterBand) = AG.getproj(AG.getdataset(a)) +# get_projection(a) = AG.getproj(a) +# get_projection(a::AbstractRasterBand) = AG.getproj(AG.getdataset(a)) function dimsymbols(a) - wkt = WellKnownText(AG.toWKT(AG.newspatialref(get_projection(a)))) - epsg = convert(EPSG, wkt) - if epsg.val == 4326 - (:Y, :X) - else - (:X, :Y) - end + # I had prepared the following code to swap x and y fpr EPSG:4326, + # However it turns out that, at least in the example https://github.com/yeesian/ArchGDALDatasets/raw/307f8f0e584a39a050c042849004e6a2bd674f99/gdalworkshop/world.tif + # dimensions are not swapped although it indicates that projection. So for now we just return x,y until we find more examples that are broken + # wkt = WellKnownText(AG.toWKT(AG.newspatialref(get_projection(a)))) + # epsg = convert(EPSG, wkt) + # @show epsg.val + # if epsg.val == (4326,) + # (:Y, :X) + # else + # (:X, :Y) + # end + (:x, :Y) end iscontdim(a::RasterDataset, i) = i < 3 ? true : nraster(a)<8 From 5d984247fa00f227b053eb004ab77d97dbef0e44 Mon Sep 17 00:00:00 2001 From: Fabian Gans Date: Thu, 21 May 2026 16:37:14 +0200 Subject: [PATCH 6/6] typo --- ext/ArchGDALExt/ArchGDALExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/ArchGDALExt/ArchGDALExt.jl b/ext/ArchGDALExt/ArchGDALExt.jl index 8f770d7..f9e29a3 100644 --- a/ext/ArchGDALExt/ArchGDALExt.jl +++ b/ext/ArchGDALExt/ArchGDALExt.jl @@ -59,7 +59,7 @@ function dimsymbols(a) # else # (:X, :Y) # end - (:x, :Y) + (:X, :Y) end iscontdim(a::RasterDataset, i) = i < 3 ? true : nraster(a)<8