Strange difference from crsTransform of subsettings from EPGS:32632 to EPSG:4326

[Bang Pham Huu]

Attachment is 2 Landsat 5 (original one is in EPSG:32632, test.tiff is warped to EPSG:4326 with gdalwarp -s_srs EPSG:32632 -t_srs EPSG:4326).

These 2 files are imported as 2 separate coverages (test_landsat5_32632 and test_landsat5_4326). If select whole coverage s in EPSG:4326, they return same result:

+ test_landsat5_4326

http://localhost:8080/rasdaman/ows?&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=test_landsat5_4326&SUBSET=ansi(%222012-12-01T20:07:00.500Z%22)&FORMAT=image/png

+ test_landsat5_32632 (needs to use crsTransform to reproject)

for c in (test_landsat5_32632) return encode(
crsTransform (
  c[ansi:"CRS:1"(0)],
{ E:"http://localhost:8080/def/crs/EPSG/0/4326", N:"http://localhost:8080/def/crs/EPSG/0/4326"}, {}
),
"png")

The results are the same (an image with 8213 x 2654 pixels).

However, if subset from these 2 coverages (coordinates from EPSG:4326 are translated to EPSG:32632 accordingly).

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32632
(Long Lat -> E N)
10 62
552375.799656895 6874583.72713382 0
11 63
601293.020582477 6987164.64881585 0

They return significantly different results.

+ test_landsat5_4326

http://localhost:8080/rasdaman/ows?&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=test_landsat5_4326&SUBSET=Lat(62,63)&SUBSET=Long(10,11)&SUBSET=ansi(%222012-12-01T20:07:00.500Z%22)&FORMAT=image/png

rasql: SELECT encode(c[1149:3168,260:2280,0], "image/png" , "{\"geoReference\":{\"crs\":\"EPSG:4326\",\"bbox\":{\"xmin\":9.999904782991448,\"ymin\":61.999546756504088,\"xmax\":11.000317768762488,\"ymax\":63.00045499622848}}}") FROM test_landsat5_4326 AS c


+ test_landsat5_32632 (needs crsTransform)

for c in (test_landsat5_32632) return encode(
crsTransform (
  c[ansi:"CRS:1"(0), E(552375.799656895:601293.020582477), N(6874583.72713382:6987164.64881585)],
{ E:"http://localhost:8080/def/crs/EPSG/0/4326", N:"http://localhost:8080/def/crs/EPSG/0/4326"}, {}
),
"png")

rasql: SELECT encode(project( c[990:2621,430:4182,0], "552354.3125,6874576,601314.3125,6987166", "EPSG:32632", "EPSG:4326" ), "png" , "{\"geoReference\":{\"crs\":\"EPSG:4326\",\"bbox\":{\"xmin\":9.999587551972840770986294955946505069732666015625,\"ymin\":61.990117737507574702249257825314998626708984375,\"xmax\":11.0005620532836463354442457784898579120635986328125,\"ymax\":63.01037371376890661167635698802769184112548828125}}}") FROM test_landsat5_32632 AS c

The first result returns an image with size 2020 x 2021 pixels. The second result returns an image with size 3011 x 3069 pixels.

Any idea about the cause?

[Bang Pham Huu]

Input files and ingredient file

[Dimitar Misev]

The subsets are translated to different grid bounds somehow:

1. First query: c[1149:3168,260:2280,0]
2. Second query: c[990:2621,430:4182,0]

It needs to be determined why this happens.

[Dimitar Misev]

Any news on comment:1?

[Bang Pham Huu]

Replying to dmisev:

Any news on comment:1?

Because they are 2 different collections with different sdoms (the first query is on an EPSG:4326 collection. The second one is on an EPSG:32632 collection).

[Bang Pham Huu]

The result is also different with using gdalwarp, so it must be due to the subsets in EPSG:4326 [Lat(62:63) and Long(10,11)] which cannot be translated point by point as [E(552375.799656895,601293.020582477 ) and N(6874583.72713382,6987164.64881585)].

Test 1:
# Subset on 32632 tiff by 32632 subset which are translated point by point from 4326 subset
gdalwarp -te 552375.799656895 6874583.72713382 601293.020582477 6987164.64881585 LT51980161984098ESA00_B1.tiff subset.tiff
Creating output file that is 1631P x 3753L.

# Warp cropped 32632 tiff file to 4326
gdalwarp -t_srs EPSG:4326 subset.tiff subset_4326.tiff
Creating output file that is 3009P x 3070L.

Test 2:
# Subset 4326 on 4326 tiff file directly
gdalwarp -te 10 61.9999999999999 11 62.9999999999954 test.tiff subset_4326_org.tif
Creating output file that is 2019P x 2019L.

[Bang Pham Huu]

I got an answer from gdal dev mailing list and opened results in QGIS as suggested

note that an EPSG:4326 rectangle of full degrees gets rotated and bended 
in almost any projected CRS (except Mercator).

Since the EPSG:32632 raster has to be an unrotated rectangle again, the 
corners of the projected EPSG:4326 rectangle are different from the 
corners of the EPSG:32632 rectangle. The difference will be filled with 
NODATA pixels.

Try to vsualize your problem in a GIS software like QGIS, and you see 
what I mean.

Basically, output in 4326 which is warped from 32632 will have more no data values in top and bottom, the warped content from 32632 also rotated while 4326 content is not.

[Dimitar Misev]

So basically reprojecting 4326 → 32632 → 4326 is lossy? We get back a result different from the original?

[Bang Pham Huu]

Replying to dmisev:

So basically reprojecting 4326 → 32632 → 4326 is lossy? We get back a result different from the original?

Yeah, it is not bi-directional possible.