SLC Incorrect Data on Write

[joshday]

I have a Single Look Complex .tiff file where the underlying data changes on a read-write-read roundtrip.

Data found on https://search.asf.alaska.edu/ (S1A_IW_SLC__1SDV_20230103T230613_20230103T230641_046626_059699_9E8C)

julia> slc = Raster(path);

julia> Rasters.write("test.tiff", slc; force=true)
"test.tiff"

julia> slc2 = Raster(ans);

# data from the "read" matrix and the "read-write-read" matrix are different
julia> slc.data == slc2.data
false

# The errors sum to zero.  An incorrect transpose?
julia> sum(x - y for (x,y) in zip(slc.data, slc2.data))
0 + 0im

julia> mean(abs(x - y) for (x,y) in zip(slc.data, slc2.data))
73.42267370323304

A likely related issue is that the written file is about half the size as the original:

julia> filesize(path)
1186109408

julia> filesize("test.tiff")
673542244

---

I'll look into this further and add more details as I find them. Just wanted to get an issue in writing.

[rafaqz]

Sounds BAD. Will try to fix. But can you download the file inline in your MWE?

Otherwise a big chunk of my available bugfix time dissapears looking for demo rasters on various platforms, rather than fixing bugs ;)

[joshday]

I know it's really annoying to reproduce...I don't think they have an API though. They require an account and all the files are > 1Gb.

My hunch is that is has to do with rotations, but I'll try to get a simpler MWE tomorrow

[rafaqz]

I just couldnt find the file after 5 mins looking a gave up 😂

[joshday]

Okay, so looks like an axis gets flipped. Here's the original and roundtrip data heatmaps side by side:

[rafaqz]

I found a file matching that name on there but it won't even load in GDAL so I'm not sure what you're actually doing - I really need a full MWE here

[joshday]

Okay, here's some better info.

@rafaqz I sent you a link to download the data. I don't think any SLC data exists in a smaller file, so that's annoying. I can't give you a smaller file but here's some simplified code that identifies the issue:

using Rasters, ArchGDAL

path = "s1a-iw1-slc-vh-20221206t234001-20221206t234026-046218-0588c1-001.tiff"

a = Raster(path)

Rasters.write("temp.tiff", a; force=true)

b = Raster("temp.tiff")  # roundtrip

a.data == b.data  # false

a.data == reverse(b.data, dims=2)  # true

[rafaqz]

Ok, this is kind of minor problem really. By convention non rotated geotiffs are usually stored with the Y axis flipped, or "north up". So users don't have to do this reversing manually all the time we do it automatically. I guess if there is any kind of rotation we should'nt do that, instead just leaving things how they are.

This is pretty well explained in this stackexchange
https://gis.stackexchange.com/questions/122647/do-rasters-other-than-north-up-exist-in-the-wild

(Thanks for the drive link its super fast to check these things with a file and MWE)