Skip to content

Only revert xy for epsg 4326 #60

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
meggart merged 6 commits into from
May 21, 2026
+68 −46
Merged

Only revert xy for epsg 4326 #60

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
64c0a48
Only revert xy for epsg 4326
meggart May 21, 2026
feb6fd5
call getproj only on datasets
meggart May 21, 2026
b77ed0c
fix some tests
meggart May 21, 2026
9b1cabc
one more try for fix
meggart May 21, 2026
bc6fb4b
remove special case again
meggart May 21, 2026
5d98424
typo
meggart May 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
102 changes: 62 additions & 40 deletions ext/ArchGDALExt/ArchGDALExt.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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)
Expand All @@ -41,6 +45,23 @@ module ArchGDALExt
end
end

# get_projection(a) = AG.getproj(a)
# get_projection(a::AbstractRasterBand) = AG.getproj(AG.getdataset(a))
function dimsymbols(a)
# 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
function getattributes(a::RasterDataset)
globatts = Dict{String,Any}(
Expand All @@ -55,45 +76,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
12 changes: 6 additions & 6 deletions test/arrays.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down
Loading