Flipped axis when plotting geographic data using `shade()`

[HTenkanen]

Hi!

Thanks for a very useful library! I am aware that there have been issues and discussion about the flipped axes when plotting xarray.DataArrays (e.g. here, here and here) and although there has been work to solve this, I believe there is still an issue when using shade() function, i.e. the output image is flipped when using the shade() function as shown below. The latter one is correct and the first one seems to be flipped:

Setup

Python 3.10.13
xarray 2023.10.1 pyhd8ed1ab_0 conda-forge
xarray-spatial 0.3.7 pyhd8ed1ab_0 conda-forge
datashader 0.16.0 pyhd8ed1ab_0 conda-forge

Is there any workaround for this, or ongoing efforts to fix this issue?

[jbednar]

There are different conventions for how to map an array onto the screen, and I think matplotlib chooses a convention that is the opposite from how arrays are normally printed, while Datashader renders in the same orientation as the data would be printed as an array. You can write xarray code to flip the data into whatever orientation you like before calling shade().

[jbednar]

(The previous issues were about having the data flipped relative to the coordinates, or flipped in some cases rather than others, which are bugs; flipping relative to some other tool's convention is not a bug. Of course, there could be a bug here as well, but I don't think the above plots demonstrate one, only a difference.)

[HTenkanen]

Hi, and thanks for a quick reply!

Hmm, I'm sorry but I do believe this is a bug (or at least something you should consider when working with geographic xarray Datasets). The data that is plotted represents an elevation map (based on a aerial/satellite image). As the image represents part of an Earth and it gets flipped, I believe it is a real problem and not only a matter of convention. In fact, DataShader used to work correctly still a couple of years ago. Meaning that this is a change in the logic of DataShader and not only a difference against "other conventions". These two maps were produced with the same shade() function but the first one (correct in geographical terms) was done with datashader=0.13 vs the latter one with current version:

I also believe I know the reason for this behavior, as this is a common mistake with my GeoIT students at University. You most likely treat Latitude and Longitude coordinates as X and Y in DataShader? If this is the case, then the axes get flipped because, in fact, Latitude=y and Longitude=x.

In xarray, which is widely used in geospatial domain the data is often repreresented with Latitude and Longitude geographic coordinates and this is the reason why x and y coordinates are "flipped". This is a peculiarity (I know) which happens because in spoken language we are very customed to talk about "Lat Lon" values, but in math those are "y and x".

I spotted this change in behavior of datashader already last year but didn't react then because I assumed someone else would spot the same issue and fix would come. Hence, somewhat surprisingly found this issue existing still and wanted to make you aware of this.

[HTenkanen]

This explanation above is the very reason for why xarray uses x-axis as dimension 1 and y-axis as dimension 0 as discussed here: holoviz/hvplot#668 (comment)

I know that not all the data is geospatial by nature, but a lot of data that is processed with xarray is. There is this mention in their documentation about this:

Hence, when plotting an image based on geographic coordinates, you should probably use the xc and yc attributes instead of x and y.

[jbednar]

To be clear, neither of the above plots involve swapping x and y; they only show inverting y (does y increase in the vertical direction, or does it decrease?). Swapping x and y results in a rotation and flip, rather than the simple vertical mirror flip as shown in these examples.

Cartesian coordinates (and geographic latitudes) have y increasing in the up direction, while computer graphics screens (and text printing of arrays) have y increasing in the down direction. Both choices are purely a matter of convention. Different xarray data files also declare their coordinates in increasing or decreasing order, which again is purely a matter of convention. What to call the coordinates is also purely a matter of convention. Datashader can't necessarily match the conventions in any particular domain or application, as it is a general rendering tool; best we can do is document what it does and how to achieve the results you want for data that follows a certain convention.

Maybe we should extend shade() to accept a name for the x and y coordinates to use, for a user who knows they want "xc" rather than "x"? shade() hasn't previously needed that because it's usually used with the output from a Datashader aggregation pipeline, which has predictable coordinate ordering and naming. I think that's a very reasonable feature request, though I'm not actually sure how to implement it because I can't see anything explicitly referencing "x" or "y" in the current shading code.

Again, I'm not saying there is no bug here, just that I haven't seen anything that demonstrates one. Differing conventions are not a bug, and the fact that it changed is also not itself a bug, because there were bugs in not respecting the xarray coordinate declarations prior to Datashader 0.15.1. Fixing those bugs will necessarily change the output, which previously did not respect what was declared in xarray.

[HTenkanen]

Thanks for the in-depth explanation! 👍🏻 Indeed, you're right, the swap in orientation between the older vs current version of datashader looks like only one axis is flipped (which indeed might be due to the earlier issues which you mentioned). I will double check the data again in another software and try to get a better understanding what is going on here.. It might be just me confused here 😛 But still my feeling is that the geographic coordinates are flipped (as in the first pair of comparison). If this is the case, then indeed being able to specify the x and y would perhaps be a nice way to deal with this. I will check the codebase to also better understand the internals.

[jbednar]

Sounds good! I believe that Matplotlib uses the Cartesian convention for plotting Numpy arrays (effectively y-inverting arrays between their printed representation and their plotted representation), but I don't know what convention xarray's .plot() uses when it calls Matplotlib (it may or may not invert it in the call), nor how it handles the various different internal orderings possible in an xarray DataArray. The whole situation of competing conventions is very confusing and we totally could have a bug in how we do it, but I'm not even sure how to figure that out if so.

[HTenkanen]

Hi @jbednar, well it took quite some time to get back to this.. But now this came timely again as I'm teaching these things this time of the year. I now finally double checked that also other GIS software (QGIS in this case) produces a hillside that looks like below:

This corresponds to the map that still worked when using datashader==0.13 (see the first comment on this issue). Currently, the orientation on the axis 0 is flipped.

I was able to overcome this issue by following/modifying the approach introduced earlier in this issue on the matter and flipping the data on axis 0 (as you suggested as well):

stack(
    shade(np.flip(data["hillshade"], 0), cmap=['gray', 'white'], alpha=255, how='linear'),
    shade(np.flip(data["elevation"], 0),     cmap=Elevation,         alpha=128, how='linear'),
    )

I'm not sure whether this is a bug in datashader but it can be confusing for people who are working with GeoTIFF files without knowing that the orientation on one axis is flipped, so at least mentioning this possibility of flipped orientation of axis in the documentation would be very useful.