Output force tile has a slightly wrong extend

[charalamm]

Hello,

I am trying to create a force datacube out of Sentinel2 but the extend of the output products is slightly off to what I am expecting to see.

For example the X0011_Y0028/20220114_LEVEL2_SEN2A_BOA.tif of my grid according to rasterio and QGIS has an extend of BoundingBox(left=540000.0, bottom=6119999.999999999, right=550000.0, top=6129999.999999999) instead of BoundingBox(left=540000.0, bottom=6120000.0, right=550000.0, top=6130000.0) although according to gdal (gdalinfo /mnt/nas/lambis/DataCube-S1/X0011_Y0028/S1_20190528_20190609_D_139.tif | grep 'Upper\|Lower') has a correct extend:

Upper Left  (  540000.000, 6130010.000) (  9d37'48.92"E, 55d18'53.48"N)
Lower Left  (  540000.000, 6120000.000) (  9d37'43.80"E, 55d13'29.68"N)
Upper Right (  550010.000, 6130010.000) (  9d47'16.67"E, 55d18'50.18"N)
Lower Right (  550010.000, 6120000.000) (  9d47'10.27"E, 55d13'26.39"N)

Does anybody know why is this happening?

I am using the following parameters among others:

DO_REPROJ = TRUE
DO_TILE = TRUE
TILE_SIZE = 10000
BLOCK_SIZE = 2500
RESOLUTION_SENTINEL2 = 10
ORIGIN_LON=7.8214115240990782
ORIGIN_LAT=57.8264679607105876
PROJECTION='PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25832"]]'

Also, on level 2 processing I am getting the following datacube-definition.prj file which looks good to me:

PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25832"]]
7.821412
57.826468
430000.000000
6410000.000000
10000.000000
2500.000000

Thanks!

[davidfrantz]

Hi @charalamm ,

sorry for the delayed answer. FORCE uses double arithmetic to compute coordinates. What you see may be related to the floating point precision or the way that QGIS handles these data types. Please see here for some explanation about it:

https://www.php.net/manual/en/language.types.float.php

In any case, I don't see any practical drawbacks of this behaviour since the difference in coordinates is extremely small (if there is any, probably just the way how QGIS prints the numbers).

Best,
David