Petascope_Different result in total pixel sizes with GDAL

[Bang Pham Huu]

Using a test coverage (in attchment) and WCPS query

for R in (test_pixel)return encode( 
(( R)[E(377983:390000),N(4902991:4917275)] ), "tiff")

it will returns a 2D trimmed image with total pixel sizes (402, 477)

but using GDALWarp to trim the same coverage

gdalwarp -te 377983 4902991 390000 4917275 
'/home/rasdaman/Downloads/ows(179)' /home/rasdaman/Downloads/output.tiff

it only has Size is 401, 476

See the diff between both results's metadata here.

​https://www.diffchecker.com/wrvtylaf

[Dimitar Misev]

Referring to #1324 I think it's actually an issue in petascope as the subset translates to [0:401,0:476] which is 402x477.

select encode(R[0:401,0:476], "GTiff",
"xmin=377983;xmax=390000;ymin=4902991;ymax=4917275;crs=EPSG:32631")
from test_pixel AS R where oid(R)=133121

[Bang Pham Huu]

Replying to dmisev:

Referring to #1324 I think it's actually an issue in petascope as the subset translates to [0:401,0:476] which is 402x477.

select encode(R[0:401,0:476], "GTiff",
"xmin=377983;xmax=390000;ymin=4902991;ymax=4917275;crs=EPSG:32631")
from test_pixel AS R where oid(R)=133121

Yes, it can be. I suspect somewhere use + 1 when translate to CellDomain Interval.

[Bang Pham Huu]

@Vlad: it has a tricky problem here (as you know Petascope core then maybe you can answer). The problem is with something called 0.5 * pixel size as origin in Petascope.

So: coverage test_pixel is geo-referenced in ESPG:32631 (UTM), 30 m / pixel.

E: 377955 390015 (0 - 401): 402 pixels
N: 4902975 4917285 (0 - 476): 477 pixels.

So if I use Gdalwarp to crop a subset which is in first below 0.5 * 30 m ( 15 m) like this

gdalwarp -te 377955 4902975 377970 4902990 '/home/rasdaman/Downloads/input.tif' /home/rasdaman/Downloads/output6.tif

Creating output file that is 0P x 0L.

It will say it is invalid subset (0 x 0 pixel), but if I increase to 1 meter in each axis, it is 1 X 1 pixel.

gdalwarp -te 377955 4902975 377971 4902991 '/home/rasdaman/Downloads/input.tif' /home/rasdaman/Downloads/output7.tif
Creating output file that is 1P x 1L.

With Petascope, it does not concern this problem and in Rasql query from WCPS 1.0, it returns the 1 X 1 Pixel as normal.

for R in (test_pixel)return encode((( R)[E(377955:377970),N(4902975:4902990)] ), "tiff")

select encode(R[0:0,476:477], "GTiff", "xmin=377955.0;xmax=377970.0;ymin=4902975.0;ymax=4902990.0;crs=EPSG:32631") from test_pixel AS R where oid(R)=133121

which has another problem (the maximum of N axis pixel is 476 (0 - 466) not 477).

in WCPS 1.5, it has other problem as the trim interval is not valid (477:476).

SELECT encode( (  ( R ) [0:0,477:476] ) , "GTiff" , 
"xmin=377955.0;xmax=390015.0;ymin=4902975.0;ymax=4917285.0;crs=EPSG:32631") FROM test_pixel AS R

Question: Should we make same behavior as GDAL (regard to 0.5 * pixel size) ? and fixed it in WCPS 1.5 (only) ?

[Bang Pham Huu]

Relation between geo-referenced bounding box and grid bounding box

[Bang Pham Huu]

Add image for idea of GDAL
​http://rasdaman.org/attachment/ticket/1325/grid_0.5_pixel.jpeg

[Bang Pham Huu]

will need to fix in ​http://rasdaman.org/ticket/1406 as Petascope has same behavior with gdalwarp.