Dear all,
I’m Sergio and from 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.
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.
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 = [info["index_x"] for info in info_list]
idxs_y = [info["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
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)
