# Copernicus Sentiel 2 overlapping coordinates between tiles

**URL:** <https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280>\
**Category:** Data\
**Created:** [March 3, 2021, 6:02pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280 "2021-03-03T18:02:52Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![pl.marasco](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/pl.marasco/32/258_2.png) [@pl.marasco](https://discourse.pangeo.io/u/pl.marasco)\
**Post date:** [March 3, 2021, 6:02pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/1 "2021-03-03T18:02:52Z")

</div>

A little bit of Context:

ESA is distributing [Sentinel 2]( [Sentinel-2 - Missions - Sentinel Online - Sentinel (esa.int)](https://sentinel.esa.int/web/sentinel/missions/sentinel-2)) A/B at level-1C (Top Of Atmosphere) and L2A (Bottom Of Atmosphere) through tiles of fixed size. These cover an area of 100x100 km and are orthorectified-images in [MGRS/WGS84](https://en.wikipedia.org/wiki/Military_Grid_Reference_System) reference system.

For historical reasons, each tile has an overlapping area with around tiles. Unfortunately, overlaps are not fixed in dimension and vary according to the tile position within the reference grid. Within the metadata, there isn’t any information about the surrounding tiles.

The classic approach:  
As most of the time, a single tile doesn’t cover the area of interest, mosaicking of the tiles is needed. Different approaches can be taken on how overlapping values are treated. Min, max, average, first, last are the most common according to the needs but much more complicate algorithms are available for specific cases.  
In Rasterio, this can be easily obtained through the [.merge](https://github.com/mapbox/rasterio/blob/522cf4d7f45b9bf3a1d9ca2f42e6b9f417cd4e8e/rasterio/merge.py) module. To have a quick example of how ti is used, have a look [here](https://automating-gis-processes.github.io/CSC18/lessons/L6/raster-mosaic.html)

Xarray:  
From my understanding seems that there is no way to achieve the merging for one specific reason. Combine can’t be used as coordinates are not monotonic. Moreover, there isn’t an option to manage the overlapping pixels, but that would be the second step once the coordinate grid is determined. Here a discussion on this topic [Xarray combine\_by\_coords return the monotonic global index error]( [Xarray combine\_by\_coords return the monotonic global index error · Issue #4213 · pydata/xarray (github.com)](https://github.com/pydata/xarray/issues/4213))

AWS:  
For those interested in why this problem seems to be not present in the AWS, what I can say is that it seems that Sinergise, for the AWS service, has converted the files from jpeg2000 to tif, making some adjustments. Unfortunatelly I don’t have direct access to the AWS and can’t check by my self, but picking some info here and there seems that they used the average algorithm. This specific approach could be considered acceptable for the L1C but not for the L2A. But this is a completely different topic, and I would not investigate more for the moment.

Question:  
Does anyone has stumbled with this problem? How did you solve it? It could be that I’m not aware of the specific module that can solve this. I’ve searched a little bit, but nothing has popped up.  
Any idea about how I could organise a workaround in a ‘lazy dask’ approach?  
It would be feasible to think of implementation in Xarray of a system to treat this specific problem.

---

<div class="post-metadata">

**Author:** ![rabernat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/rabernat/32/22_2.png) [@rabernat](https://discourse.pangeo.io/u/rabernat)\
**Post date:** [March 4, 2021, 3:31pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/2 "2021-03-04T15:31:27Z")

</div>

Good question @pl.marasco!

@JessicaS11, @scottyhq and @TomAugspurger, and others have been working on developing some best practices for working with geospatial raster data in Pangeo. This is discussed in the following threads

> [@Xarray for raster data (DEMs) with inconsistent spatial extent](https://discourse.pangeo.io/t/xarray-for-raster-data-dems-with-inconsistent-spatial-extent/821):
>
> I’ve been spinning up a new project and have been trying to use xarray, but I’m remembering why I abandoned it years ago the first time I tried to use it. It’s not intuitive to teach yourself (e.g. a DataArray automatically has a variable included that’s not listed in the output, but a DataSet lists the variables by name, even if there’s only one), and the examples/resources are thin to nonexistent if your data don’t fit into nice, neat consistent grids (whether that’s because they have differen…

> [@Cloud Optimized Geotiffs + Pangeo best practices](https://discourse.pangeo.io/t/cloud-optimized-geotiffs-pangeo-best-practices/999):
>
> Hi everyone, I’m attempting to consolidate collective wisdom on working with Cloud Optimized Geotiffs on Pangeo infrastructure. Mainly focusing on using Xarray and Dask to efficiently run computations against collections of COGs from within the same Cloud datacenter. I think there is quite a bit to explore here, so I’ve decided to put together a repository to add several notebooks to get started. More details in the repository readme: I think there are already some helpful resources and exa…

Perhaps they will have some suggestions.

---

<div class="post-metadata">

**Author:** ![c00kiemon5ter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/c00kiemon5ter/32/578_2.png) [@c00kiemon5ter](https://discourse.pangeo.io/u/c00kiemon5ter)\
**Post date:** [March 4, 2021, 6:45pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/3 "2021-03-04T18:45:09Z")

</div>

hello, I agree with @rabernat - good question!

I think this is a super hard problem to solve, and tbh, I don’t think it is solvable at the level of the products given out by ESA.

Think of this; ESA hands out tiles - squares based on MGRS - and in that it makes a simplification. The assumption is that, given that dimension of a tile, the curvature of Earth is not affecting the area representation. What you are trying to do is to stitch together multiples of those squares, but suddenly, when you start doing this the dimensions grow, and curvature does start to play a role. The southern or northern you work, the more prevalent that is.

I don’t think there is a good way to stitch together L1C or L2A tiles on a 2D space; at least there doesn’t seem to be one because it would not result in a useful representation (or visualization even).

(I will also leave this here which is not answering your question, but it is interesting:)

> **[Changing the way you look at Earth data with xESMF:](https://medium.com/pangeo/changing-the-way-you-look-at-earth-data-with-xesmf-ee55d0b380e5)**
>
> How our different ways of looking at earth data became a mathematical and software problem

Having said this, I am assuming that you are trying to solve a different problem, and this is a step in the process to solve it. What is it that you are really trying to do?  
In most cases, you would be processing the image to understand features and deduct something about the given area. In most cases you want to avoid processing the same area (part of tile) twice, both due to performance reasons, but also due to the bias of evaluating such parts multiple times.  
Instead of thinking of a big surface, composed of square images, try to find for each pair of neighbouring tiles those common pixels. Then think of processed and unprocessed areas/pixels, guided by this information.

---

<div class="post-metadata">

**Author:** ![rabernat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/rabernat/32/22_2.png) [@rabernat](https://discourse.pangeo.io/u/rabernat)\
**Post date:** [March 4, 2021, 6:50pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/4 "2021-03-04T18:50:37Z")

</div>

> [@pl.marasco](#):
>
> In Rasterio, this can be easily obtained through the [.merge](https://github.com/mapbox/rasterio/blob/522cf4d7f45b9bf3a1d9ca2f42e6b9f417cd4e8e/rasterio/merge.py) module. To have a quick example of how ti is used, have a look [here](https://automating-gis-processes.github.io/CSC18/lessons/L6/raster-mosaic.html)

Rereading your original post, it sounds like a good solution already exists. Can you help us understand better why the rasterio solution is not enough for you?

---

<div class="post-metadata">

**Author:** ![JessicaS11](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/jessicas11/32/229_2.png) [@JessicaS11](https://discourse.pangeo.io/u/JessicaS11)\
**Post date:** [March 4, 2021, 8:02pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/5 "2021-03-04T20:02:41Z")

</div>

> [@pl.marasco](#):
>
> Xarray:  
> From my understanding seems that there is no way to achieve the merging for one specific reason.

Will the [`rioxarray.merge` module](https://corteva.github.io/rioxarray/stable/rioxarray.html#rioxarray-merge-module) (which implements `rasterio.merge`) work for your purposes, or are you looking for a way to merge the rasters that does not use rasterio?

---

<div class="post-metadata">

**Author:** ![rmg55](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/rmg55/32/971_2.png) [@rmg55](https://discourse.pangeo.io/u/rmg55)\
**Post date:** [March 4, 2021, 8:24pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/6 "2021-03-04T20:24:35Z")

</div>

@pl.marasco - in my experience, this is a common challenge when developing workflows with these data.

My current approach is to leverage [gdal vrt](https://gdal.org/drivers/raster/vrt.html) to stack (temporally) or stitch(mosaic) these data. The vrt file can also specify warp parameters ([gda\_vrt\_warping](https://gdal.org/drivers/raster/vrt.html#warped-vrt)) to reproject to a new crs (e.g. epsg 4326).

@TomAugspurger has just started an effort to make this process easier, which leverages STAC catalogs ([GitHub - stac-utils/stac-vrt](https://github.com/stac-utils/stac-vrt)). I think at the moment it focuses on the mosaic workflow (example with NAIP imagery), but I assume a stacking (temporal) functionality can be implemented.

---

<div class="post-metadata">

**Author:** ![pl.marasco](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/pl.marasco/32/258_2.png) [@pl.marasco](https://discourse.pangeo.io/u/pl.marasco)\
**Post date:** [March 5, 2021, 8:37am UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/7 "2021-03-05T08:37:45Z")

</div>

Hi @c00kiemon5ter,

Projection issue:  
Tiles are projected in a 2d dimension through a Tranverse of Mercator to be stitchable within a zone. Simplifying quite a lot this type of projection, conformal ones, preserves the ratio of two lengths (in this case angle and length ) keeping errors acceptable within a ‘small’ spatial domain. You would not have a problem to create a single huge image if you will stitch within the boundary of 6 longitude degree.  
Indeed, once some on will try to deal with two different zones, a completely different approach would be taken but I would not consider it now.

xESMF:  
I had a quick look at it and seems really interesting!!! thanks a lot, indeed I’m going to use it. The best way to treat different zones is to convert from planner coordinate to geographic once and the xESMF seems to treat exactly these issues. Unfortunatelly the overlapping is not going to be solved even with this change of reference system.

What you described is correct. The only problem is that once you will start to create, in practical therm, such thinks the first problem that you will figure out is that Xarray can’t deal with overlaps. Consequentially you will have to find a smart way to deal with it keeping in mind that without a lazy approach will end up with a huge amount of data.

Target:  
Create a cube of data where I’ve, for each sensed band (let’s say Red/Green/Blu) a mosaic over wich I can apply some algorithms that will involve the band and the time dimension.

---

<div class="post-metadata">

**Author:** ![pl.marasco](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/pl.marasco/32/258_2.png) [@pl.marasco](https://discourse.pangeo.io/u/pl.marasco)\
**Post date:** [March 5, 2021, 8:53am UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/8 "2021-03-05T08:53:55Z")

</div>

@rabernat  
Indeed .merge is the usual solution for this type of things but I was trying to create a workflow based on a completely lazy approach so the final dataset will manageable.

In any case, I found it interesting that the combine by coords can’t deal with this and I thought to bee a good idea to raise the issue in the community. More easily users can deal with their data more users will find it convenient to use it. Real few people rise the issue as most of the others are leveraging on the work of Synergise and AWS.

Even if we are in a testing phase the entire area covered by the project is bounded within West Africa to India, temporal resolution speaking is 120 days. So I have to create something that will be able to deal with a huge amount of data and I was more prone to create something that not duplicate data on the disk (or at least not from the first step)

@JessicaS11

Good point! I wasn’t aware of rioxarry and seems that most probably will solve the problem! I’ll let you all know if it was fitting the purpose.

@rmg55  
Indeed VRT is a good solution to create virtual raster and I’m used too it. Unfortunatelly the wrapping would require a huge amount of effort as you would need for each tile know the neighbour tiles and the exact non-overlapping area. That’s why on a different group the first question has been … ‘Did ESA provide, within the metadata, the non-overlapping area?’ unfortunately I’ll let you imagine the answer.

---

<div class="post-metadata">

**Author:** ![bluetyson](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/bluetyson/32/763_2.png) [@bluetyson](https://discourse.pangeo.io/u/bluetyson)\
**Post date:** [March 6, 2021, 10:01pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/9 "2021-03-06T22:01:43Z")

</div>

Vrt has the problem of dealing with multiple crs…so if a large area, you have to handle that first.

---

<div class="post-metadata">

**Author:** ![RichardScottOZ](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/richardscottoz/32/752_2.png) [@RichardScottOZ](https://discourse.pangeo.io/u/RichardScottOZ)\
**Post date:** [March 6, 2021, 10:06pm UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/10 "2021-03-06T22:06:37Z")

</div>

The simple merges in rasterio may be insufficient for spectral reasons.

---

<div class="post-metadata">

**Author:** ![pl.marasco](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/pl.marasco/32/258_2.png) [@pl.marasco](https://discourse.pangeo.io/u/pl.marasco)\
**Post date:** [March 8, 2021, 6:54am UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/11 "2021-03-08T06:54:34Z")

</div>

Good point, if someone needs to cover a large area that is span over more than one MRGS zone one of the best options would be to convert everything in lat/lon and then merge.  
Indeed would be nice to know how Synergise/AWS treated the projection. Did they move to a Geographic one? From the FAQ seems that for a bug, the L2A lost the metadata that brings the CRS info (no idea if they fix it meanwhile).

---

<div class="post-metadata">

**Author:** ![pl.marasco](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/pl.marasco/32/258_2.png) [@pl.marasco](https://discourse.pangeo.io/u/pl.marasco)\
**Post date:** [March 8, 2021, 7:24am UTC](https://discourse.pangeo.io/t/copernicus-sentiel-2-overlapping-coordinates-between-tiles/1280/12 "2021-03-08T07:24:14Z")

</div>

I wouldn’t be so strict in the answer. You would need to know what you are doing according to the specific project.  
Indeed would be a mistake use an avarage on a [L2A](https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/product-types/level-2a) scene; at the same time I would not be so restrictive using a First or Last approach in that specific case.  
For [L1C](https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/product-types/level-1c) the situation should be solved using an algorithm that fits best the purpose (min/max/avarage…)

Thanks to the really interesting and well done [notebook](https://github.com/RichardScottOZ/dask-era5/blob/main/notebook/PANGEO-Sentinel-Intake-SA-OZLocation.ipynb) that you shared few weeks ago my attention has been focus on the metadata from the AWS products. ‘OVR\_RESAMPLING\_ALG: AVERAGE,’ is reported in a L2A product. I’ve been searching for a little in the documentation but I’ve not been able to fid anything regarding this. Maybe you have more info on this.  
For sure, if that would be the algorithm, I would not be too happy as L2A is a [BRDF](https://en.wikipedia.org/wiki/Bidirectional_reflectance_distribution_function) correct and in my application would be better to have original values (First/Last alg.).