TileSplitter function problems with Sentinel2 images

Dear all,
My company I’m training to use sentinelhub services and api.

Specifically, I am using your API to locate and download Sentinel2 image projections of interest to me using the eolearn library.

I’ve encountered a problem in my workflow and wanted to try a question. The scenario is as follows:

I locate an area of interest
I create the grid of the patch of interest with the API “TileSplitter” correctly setting the CRS to EPSG:4326 (which in addition to being the crs of sentinel2 is also the crs of the polygon that describes my area of interest).

As you can see from the image the area in red is my area of interest, while in green are shown the output patch from the function TileSplitter.

finer_tile_splitter  = TileSplitter(
    rboundary_shape],
    CRS.WGS84,
    ('2017-11-1', '2020-11-2'),
    data_collection= DataCollection.SENTINEL2_L2A,
    config=config,
    tile_split_shape =6,
)

From this operation Note that some tiles are neither adjacent nor overlapping and therefore leaving areas not covered by any image and this for my type of application is a problem.

Am I doing something wrong ? Or is it an API problem ? If it is an API problem can I fix it?

Thank you very much

Page 1 / 1

Hi,

Bounding boxes of Sentinel-2 tiling grid are originally in UTM coordinate reference systems. The “correct” transformation of such bounding boxes into WGS 84 are not bounding boxes anymore, they are skewed and slightly rotated shapes. Therefore, any implementation that transforms a bounding box from one CRS into a bounding box in another CRS is only an approximation. So the inaccuracies that you see when you plot bounding boxes in WGS84 are because somewhere in your process such an approximate transformation is used.

The place where this happens might be if you call:

tile_splitter.get_bbox_list(crs=CRS.WGS84)

Instead I suggest you call:

tile_splitter.get_bbox_list()

This will provide bounding boxes in original UTM CRSs and you should continue working with those directly.

Thanks for your answer sir. So, Is not possibile to convert BBox in EPSG:4326 without having this type of issues ?

config = SHConfig()

if config.instance_id == '':
    print("Warning! To use WFS functionality, please configure the `instance_id`.")

finer_tile_splitter  = TileSplitter(
    shape_list==boundary_shape],
    crs=getCRS.from_epsg("4326"),
    time_interval=('2017-11-1', '2020-11-2'),
    data_collection= DataCollection.SENTINEL2_L2A,
    config=config,
    tile_split_shape =6,
)

bbox_list = np.array(finer_tile_splitter.get_bbox_list(getCRS.from_epsg('4326')))
info_list = np.array(finer_tile_splitter.get_info_list())
geo =  geometry.Polygon(bbox.get_polygon()) for bbox in bbox_list]
idxs_x =  infoo"index_x"] for info in info_list]
idxs_y =  infoo"index_y"] for info in info_list]
tile_id =  from_splitInfoList_to_listOfTilesID(info) for info in info_list]

bbox_gdf = gpd.GeoDataFrame({"tile_id":tile_id, "index_x": idxs_x, "index_y": idxs_y}, crs=boundary.crs, geometry=geo)

# Creation of shp file that show selected tiles
shapefile_name = "grid_aoi.shp"
bbox_gdf.set_crs("EPSG:4326")
bbox_gdf.to_crs("EPSG:4326")
bbox_gdf.to_file(os.path.join(RESULTS_FOLDER, shapefile_name))

Is i understood correctly your answer the problem is in

bbox_list = np.array(finer_tile_splitter.get_bbox_list(getCRS.from_epsg('4326')))

and your solution is to use:

bbox_list = np.array(finer_tile_splitter.get_bbox_list())

but in this way the bbox is in EPSG:32632 projection and i’m not able to translate in EPSG:4326 without distrosion.

You can if you transform them as geometries/polygons instead bounding boxes. Using helper utilities from sentinelhub-py you can do it like this:

from sentinelhub import Geometry, CRS



bboxes = finer_tile_splitter.get_bbox_list()

utm_geometries = [Geometry(bbox, crs=bbox.crs) for bbox in bboxes]

wgs84_geometries = [geo.transform(CRS.WGS84) for geo in utm_geometries]

shapely_geometries = [geo.geometry for geo in wgs84_geometries]



bbox_gdf = gpd.GeoDataFrame(..., geometry=shapely_geometries)

bboxes = finer_tile_splitter.get_bbox_list()
utm_geometries =  Geometry(bbox, crs=bbox.crs) for bbox in bboxes]
wgs84_geometries =  geo.transform(CRS.WGS84) for geo in utm_geometries]
shapely_geometries =  geo.geometry for geo in wgs84_geometries]

Thank you very much sir !