ArcticSnow / TopoPyScale

TopoPyScale: a Python library to perform simplistic climate downscaling at the hillslope scale
https://topopyscale.readthedocs.io
MIT License
39 stars 9 forks source link

handling projection of ERA5 data from lat/lon to other projection #66

Closed ArcticSnow closed 1 year ago

ArcticSnow commented 1 year ago

Currently the method to convert ERA5 grid cell coordinates from one projection system to another one is flawed. It is OK for small scale project but for larger scale we need to perform an interpolation for reprojecting.

https://github.com/ArcticSnow/TopoPyScale/blob/aa27ff1ce331bbac6fde20e4cecadee95b077298/TopoPyScale/topo_scale.py#L159

Suggested method to overcome the problem:

  1. create a grid in the target CS
  2. convert this grid to lat/lon CS
  3. sample ERA5 dataset `ds.sel(lon=xs, lat=ys, method='interp')
  4. save result to new dataset.

alternative: work in lat/lon until horizontal interpolation.

https://github.com/ArcticSnow/TopoPyScale/blob/aa27ff1ce331bbac6fde20e4cecadee95b077298/TopoPyScale/topo_scale.py#L186

replace row.x, row.y to row.lon, row.lat

ArcticSnow commented 1 year ago

Commit https://github.com/ArcticSnow/TopoPyScale/commit/0bb2e2489600913306fb749a898f413316801bff should provide a solution. @joelfiddes can you check on your side it makes sense for large scale project like yours? I used the alternative solution.

krisaalstad commented 1 year ago

It's certainly worth looking more into how we combine the projections and interpolation. I haven't thought about it in a while, but to me it always made the most sense to do spatial interpolation in a projected coordinate system or at least use something like great-circle distance if you work in regular geographic coordinates. I see that (e.g.) UTM projections can be an issue when working over very large scales, so it's certainly not a panacea. Still, I would just caution about interpolating ERA5 based on latitude and longitude directly, althought that's maybe not where this was going.

ArcticSnow commented 1 year ago

So until now TopoPyScale directly converted the grid coordinates of era5 to projected coordinate system (CS). The way it was done was wrong, not allowing for wrapping (interpolate) the grid to the new projection. With the fix done yesterday, coordinates are converted only for the point at which downscaling occurs. That way we do not need to handle projection of ERA5 at all. I made sure distance between ERA5 gridcell centroids are computed in a projected CS. So no issue there as I understand it. To conclude, the DEM at the moment must be in a projected CS. For very large scale project, then it will require to split in sub projects, or implement the method you're pointing out @krisaalstad .