ks905383 / xagg

Aggregating gridded data (xarray) to polygons
https://xagg.readthedocs.io/
GNU General Public License v3.0
83 stars 15 forks source link

Longitude -180 pixels included in weightmap #51

Closed jrising closed 9 months ago

jrising commented 1 year ago

I am aggregating ADM1-level regions (states) across the globe according to the ERA5 map. I discovered that my temperatures where, in some cases, much lower than they should be. This is because pixels from the international dateline are being included in their aggregations. I have a partial diagnosis for this problem. By the way, this is all under xagg v0.3.0.2, since I get an error (ValueError: new dimensions ('lat', 'bnds') must be a superset of existing dimensions ('lat', 'bnd')) when I try to run xagg v.0.3.1.

For example, Oslo (as defined in GADM 3.6) just overlaps 4 ERA5 pixels. However, the weight map has significant weights on 6 pixels: [(59.75, -180.0), (60.0, -180.0), (59.75, 10.75), (60.0, 10.5), (60.0, 10.75), (60.0, 11.0)]. The first two of these should not be there.

The pixel polygon that overlaps Oslo and produced the first pixel entry has these coordinates: [(179.875, 59.625), (179.875, 59.875), (-179.875, 59.875), (-179.875, 59.625), (179.875, 59.625)] which in turn comes from the first ERA5 longitude value having bounds of [ 179.875, -179.875].

So, geopandas interprets that as a strip across the whole world, which is then assigned the temperature at 180 W.

I would expect this to be a problem everywhere, but (for reasons I haven't figured out), it seems to mainly be a problem in Norway and China.

ks905383 commented 1 year ago

Yikes, that should've theoretically been captured by test_get_bnds_wraparound(). Also really weird that it's not all pixels with ERA5, given that it looks like the issue is datasets with pixel locations directly at 180/-180... Let's see what we can do.

ks905383 commented 10 months ago

Okay I think I've isolated the issue. It happens when pixels cross the antimeridian and is related to the fact that shapely doesn't know about the crs and therefore doesn't realize that coordinates should roll over at 180/-180.

The fix seems to be to change the way pixel polygons are made at the antimeridian, maybe akin to what the antimeridian package does - ie, split pixels at the antimeridian and recombine them as a multipolygon instead. Hopefully can fix soon.