Open-EO / openeo-geopyspark-driver

OpenEO driver for GeoPySpark (Geotrellis)
Apache License 2.0
26 stars 4 forks source link

OpenEO produced image with spatial distortion and discontinuities as compared with original image #829

Open sebbelese opened 3 months ago

sebbelese commented 3 months ago

Hello,

I downloaded a subset of a specific Sentinel-2 product with openEO, and also extracted the same subset directly from the .SAFE product from Copernicus dataspace.

It Generally looks very similar, and the values match, e.g. true color below.

However, if we zoom, we see that there are changes in the positioning of the features. OpenEO seem to show some discontinuities, suggesting that the issue is on OpenEO processing.

From SAFE product FromOpenEO

I hope that the difference is visible enough. It can clearly be seen when switching layers in e.g. QGIS.

Here is a minimal code to download the openEO tile. The images are also attached (note that openEO image has only 3 bands while I kept all bands in the other one from the SAFE product). The same issue occur with all sentinel images that I downloaded.

import openeo
from shapely import wkt
import datetime
from openeo.processes import eq

connection = openeo.connect("https://openeo.dataspace.copernicus.eu/openeo/1.2")
connection.authenticate_oidc()

TILE_ID = "04QHG"
CRS = "EPSG:32605"
S2_PRODUCT_ID = "S2A_MSIL2A_20190116T210921_N0211_R057_T05QKB_20190116T223213"
poly = wkt.loads("POLYGON ((222380 2123230, 227480 2123230, 227480 2128330, 222380 2128330, 222380 2123230))")
bounds = poly.bounds

sensing_time = S2_PRODUCT_ID.split("_")[2]
date_start = datetime.datetime.strptime(sensing_time,'%Y%m%dT%H%M%S')
date_end = date_start + datetime.timedelta(seconds=1)
cube = connection.load_collection(
    "SENTINEL2_L1C",
    temporal_extent=[date_start, date_end],
    bands= ['B02', 'B03', 'B04'],
    spatial_extent = {
        "west": bounds[0],
        "south": bounds[1],
        "east": bounds[2],
        "north": bounds[3],
        "crs" : CRS
    },
    properties = {"tileId": lambda x: eq(x,TILE_ID)}
)

DN_reflectance_factor = 0.0001
cube = cube * DN_reflectance_factor

result = cube.save_result("GTiff")
job = result.create_job()
job.start_and_wait(lambda a : None)
job.get_results().download_file("test_openEO.tif")

images.zip

jdries commented 3 months ago

I indeed see some things which look like different warping, so we need to check some details:

Can you pick a utm zone/tile id, and make sure to use it consistently? This can make a difference for these very small artifacts. In the result metadata of openEO, you should also see a 'canonical link'. That allows you to share the result with metadata here, and gives us more insight in which source products was used in openEO. If you can share this info, we might be able to tell where the difference is coming from.

sebbelese commented 3 months ago

Thanks for your answer. Indeed, if the tile is in the same CRS as the one used to determine the bounds (as in the code below), the alignment is correct.

These product names, tile ids and bounds come from an external database of labelled images to be used for training, and I am not sure why there are these inconsistencies in the CRS.

With the tile id below, I have well-aligned images.

import openeo
from shapely import wkt
import datetime
from openeo.processes import eq

connection = openeo.connect("https://openeo.dataspace.copernicus.eu/openeo/1.2")
connection.authenticate_oidc()

TILE_ID = "05QKB"
CRS = "EPSG:32605"
sensing_time = "20190116T210921"

poly = wkt.loads("POLYGON ((222380 2123230, 227480 2123230, 227480 2128330, 222380 2128330, 222380 2123230))")
bounds = poly.bounds

date_start = datetime.datetime.strptime(sensing_time,'%Y%m%dT%H%M%S')
date_end = date_start + datetime.timedelta(seconds=1)
cube = connection.load_collection(
    "SENTINEL2_L1C",
    temporal_extent=[date_start, date_end],
    bands= ['B02', 'B03', 'B04'],
    spatial_extent = {
        "west": bounds[0],
        "south": bounds[1],
        "east": bounds[2],
        "north": bounds[3],
        "crs" : CRS
    },
    properties = {"tileId": lambda x: eq(x,TILE_ID)}
)

DN_reflectance_factor = 0.0001
cube = cube * DN_reflectance_factor

result = cube.save_result("GTiff")
job = result.create_job()
job.start_and_wait(lambda a : None)
job.get_results().download_file("test_openEO.tif")

The canonical link found with job.get_results().get_metadata()["links"][2]["href"] for this working trial is https://openeo.dataspace.copernicus.eu/openeo/1.2/jobs/j-240613cd92b84e5488fa8d2f4147a21c/results/Njk5MDRhOGItOTJlZi00YTc0LTg4NTMtODI2Zjc3MzhlY2Y3/161678e8e9ae3d33d76e713dad05397e?expires=1718884270

Here is the canonical link for the first run with issues https://openeo.dataspace.copernicus.eu/openeo/1.2/jobs/j-240613e551994e45ac63868be9dab1ac/results/Njk5MDRhOGItOTJlZi00YTc0LTg4NTMtODI2Zjc3MzhlY2Y3/78308338b53f3a91a1b00b9174577706?expires=1718884660