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, see: ​http://rasdaman.org/attachment/ticket/1882/display_result.png

[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.