Xarray combine_by_coords return the monotonic global index error

[hamiddashti]

I asked this question on Stackoverflow and someone mentioned it might be a bug.

I am trying to combine two spatial xarray datasets using combine_by_coords. These two datasets are two tiles next to each other. So there are overlapping coordinates. In the overlapping regions, the variable values of one of the datasets is nan.

I used the "combine_by_coords" with compat='no_conflicts' option. However, it returns the monotonic global indexes along dimension y error. It looks like it was an issue before but it was fixed (issue 3150). So I don't really know why I get this error. Here is an example (the netcdf tiles are here):

import xarray as xr

print(xr.__version__)
>>>0.15.1

ds1=xr.open_dataset('Tile1.nc')
ds2=xr.open_dataset('Tile2.nc')
ds = xr.combine_by_coords([ds1,ds2], compat='no_conflicts')
>>>...
 ValueError: Resulting object does not have monotonic global indexes along dimension y

[TomNicholas]

Hi @hamiddashti , based on your description then this isn't a bug, it's throwing the error it should be throwing given your input. However I can now see how the documentation doesn't make it very clear as to why this is happening!

combine_by_coords and combine_nested do two things: they concatenate (using xr.concat), and they merge (using xr.merge). The concatenate step is never supposed to handle partially overlapping coordinates, and the combine functions therefore have the same restriction.

That error is an explicit rejection of the input you gave it: "you gave me overlapping coordinates, I don't know how to concatenate those, so I'll reject them." Normally this makes sense - when the overlapping coordinates aren't NaNs then it's ambiguous as to which values to choose.

In your case then you are asking it to perform a well-defined operation, and the discussion in the docs about merging overlapping coordinates here implies that compat='no_conflicts' would handle this situation. Unfortunately that's only for xr.merge, not xr.concat, and so it doesn't apply for combine_by_coords either. merge deals with variables of the same size, concat joins variables of different sizes onto the ends of one another. This is definitely confusing.

It might be possible to generalise the combine functions to handle the scenario you're describing (where the overlapping parts of the coordinates are specified entirely by the non-NaN values), but I don't really think it's a good idea for complexity reasons.

(Issue #3150 was about something else, an actual bug in the handling of "coordinate dimensions which do not vary between each dataset".)

Instead, what you need to do is trim off the overlap first. That shouldn't be hard - presumably you know (or can determine) how big your overlap is, and all your NaNs are on one dataset. You just need to use the .isel() method with a slice. Once you've got rid of the overlapping NaNs then you should be able to combine it fine (and you shouldn't need to specify compat either). If you're using combine_by_coords as part of opening many files with open_mfdataset then it might be easier to write a trimming function which you apply first using the preprocess argument to open_mfdataset.

Would it be clearer if I added an example to the documentation where some overlapping data had to be trimmed first?

[hamiddashti]

Thanks @TomNicholas for a thorough explanation. Now it makes sense. I thought this is a process like mosaicing rasters using rasterio.merge. Yes, it would be great if you can add an example on how to find and trim overlapping coordinates. I don't really know how to find common coordinates in spatial datasets using isel.

Thanks

[TomNicholas]

Thanks @TomNicholas for a thorough explanation. Now it makes sense.

No problem, thanks for flagging it.

Yes, it would be great if you can add an example on how to find and trim overlapping coordinates. I don't really know how to find common coordinates in spatial datasets using isel.

Do you already know which data points overlap? You know they are NaNs, so do you know how many NaNs there are at the edge of your tiles? If you do then it's just like .isel(lat=slice(2, None)) for example, but if you don't then it might be a more complicated pre-processing to do.

[TomNicholas]

What I meant is that .isel doesn't find common coordinates, it just selects points based on index. Selecting based on coordinate is done with .sel. And if you want to interrogate the values of the coordinates first you will need to do some kind of comparison or logical selection.

[hamiddashti]

Thanks. The problem is I don't really know where the common coordinates are. However, since I know the extent of each tile, with some preprocessing I should be able to find/trim them.

On the side note. It would be great if mosaicing dataset with common coordinate was added to the xarray (like mosacing rasters). What really happened was I had a large area, then I used xarray to clip it into smaller tiles to make calculations more feasible. After I did the calculations on tiles and decided to stitch them together this problem raised. Your approach would solve this issue, but having it as part of xarray definitely helps many geospatial applications where, for many reasons, we deal with tiles of data with common coordinates.
Thanks again.

[TomNicholas]

Thanks. The problem is I don't really know where the common coordinates are. However, since I know the extent of each tile, with some preprocessing I should be able to find/trim them.

Great. Let me know if you still have problems (on here, SO - I just answered your original question there, or on the xarray mailing list).

What really happened was I had a large area, then I used xarray to clip it into smaller tiles to make calculations more feasible.

I wonder if you could have avoided having to do this by applying your analysis in chunks using dask? That might be complicated if your analysis is a complicated algorithm though.

It would be great if mosaicing dataset with common coordinate was added to the xarray (like mosacing rasters). ... having it as part of xarray definitely helps many geospatial applications where, for many reasons, we deal with tiles of data with common coordinates.

This sounds like something that might be useful for lots of geoscientists, so it would be good to discuss this further. However, I don't really know exactly what you mean by "mosaicing rasters" (I don't work in geoscience). Briefly reading about it here it seems that there isn't one universal way to do it... What would be really great is if you could give me a more precise specification of the behaviour you're imagining, and how it would be used in practice (either here on in a new issue). Then we can see if it's (a) feasible, (b) commonly-useful, and (c) should live in xarray or another package. Another good place to ask about the best way to approach this problem in general would be the pangeo discourse.

[pl-marasco]

@TomNicholas I've landed on this discussion looking for a solution for what I consider the exact same problem.
Indeed the overlapping is something that all the users of Sentinel 2 Level 1c will figure out.
All the observations are deployed to users through a series of tiles following the MGRS grid system. Each tile has an overlapping area with the bordered once and is varying according to the position of the tile in relation to the reference system.
Indeed the approach you are describing can solve the problem but would require the analysis of the bounding box and a consequential selection through the .sel().
In Rasterio this can be easily obtained through the .merge module.
To have a quick example of how ti is used have a look here

You are right in pointing that there are multiple ways to treat the overlapping values but I would stick with the most common one that is as well reported in the link you mentioned. In other words (min, max, average, first, last) would be already a huge plus.

About dask, indeed is helping a lot to create a delayed object of the tiles (consider that at least for S2 data are in jp2 and we are forced to use open_rasterio instead of open_mfdataset) so the solution should be compatible with this kind of approach.
If you need further explanation or I wasn't too clear please let me know.

About Pangeo, indeed a topic should be opened on it and eventually we can move there the discussion but, at least in my opinion, for the moment the right place to discuss is within xarray.

Seems that Sinergise for the AWS service has used the average algorithm to solve the same issue. Seems that all the users that will use the AWS S2 Products will not need to care about the overlap issue.

Edit: update on AWS service

[TomNicholas]

Hi @pl-marasco, thanks for your comment.

So what you're suggesting is to alter combine_by_coords and combine_nested to be able to optionally handle overlapping data, presumably with an additional keyword argument to specify the rule it uses?

I think that this could be done within the combine functions, but it's not trivial... At the moment if you pass data with big gaps in it to combine_by_coords it can fill it with NaNs, so that's at least sort of the inverse of what you're suggesting.

Internally the combine functions currently work by creating an intermediate representation of the arrangement of tiles, before combining that along 1D repeatedly until done. What I'm wondering is whether any treatment of overlapping values would need to happen before this 1D combining step? If I have 4 tiles which all overlap at a corner, and you want me to take the (max/min/average) value of all 4 in that quadruple overlap region, I could either do this by identifying that region and taking the max (complicated) or by simply updating the max value every time I combine along 1D (simple, but more wasteful). Separately, a treatment based on the order of the input passed (your first/last) would I think need to store extra information about that order, which would be more complicated.

Do these raster problems always use the same sized tiles?

[TomNicholas]

@pl-marasco I've just been pointed to this issue on pangeo-data, which looks like a better place to discuss this.

[TomNicholas]

Closing this as having answered the original question. If anyone wants to discuss mosaicing rasters in more detail we can raise another issue for that.