JiaweiZhuang / xESMF

Universal Regridder for Geospatial Data
http://xesmf.readthedocs.io/
MIT License
269 stars 49 forks source link

regridding yields 0 when there should be no value #63

Closed pascalgross closed 4 years ago

pascalgross commented 4 years ago

I use the following code to:

  1. read a grib file with weather data
  2. regrid it using xesmf
  3. export it as GeoTIFF using rioxarray
import cfgrib
import rioxarray
import xarray as xr
import xesmf as xe
import numpy as np
from src.datapaths import get_grib_file
from src.mkdir_p import mkdir_p

year, month, day, hour = 2015, 5, 10, 12
grib_file = get_grib_file(year=year, month=month, day=day, hour=hour)
datasets = cfgrib.open_datasets(grib_file, backend_kwargs={'indexpath': ''})
ds = datasets[0]

ds = ds.rename(latitude='lat', longitude='lon')
ds_out = xr.Dataset({'lat': (['lat'], np.linspace(44.77, 56.14, 1000)),
                     'lon': (['lon'], np.linspace(2.976, 19.84, 1000))})
regridder = xe.Regridder(ds, ds_out, 'bilinear', reuse_weights=True)
ds_lonlat = regridder(ds)
ds_lonlat = ds_lonlat.rename(lat='y', lon='x')

for var in ds_lonlat.variables:
    if ds_lonlat[var].dims == ('y', 'x'):
        ds_lonlat[var].rio.set_nodata(9999).rio.set_crs(4326).rio.to_raster('{}.tif'.format(var))

This works fast as hell and is exactly what I needed. With one problem: values which are not available in the source are set to 0.0

The results in "shadows" for values outside the source grid. This is the actual result and this is what I expected.

Can I somehow define, which value is used, if no value can be interpolated?

JiaweiZhuang commented 4 years ago

values which are not available in the source are set to 0.0

Yes, this behavior is discussed in #15 #51

In short, "regridding" is just a sparse matrix multiply; the output matrix from a sparse matrix multiply is zero by default.

Can I somehow define, which value is used, if no value can be interpolated?

Filling unmapped points with non-zeros values will make the entire regridding operation non-linear (not like y = A*x anymore), so it is not possible to enable a non-zero default by just changing the regridding weights. One hack I can think of is to use https://github.com/JiaweiZhuang/xESMF/issues/15#issuecomment-371646763 to convert unmapped points to nan and then convert nan to a user-specified value. Or, use masks (#22) instead of hacking nans. Welcome more clever suggestions...

A side question: how did you make the map plots? They look quite fancy...

pascalgross commented 4 years ago

Okay, thank you. I'll have a look on that.

The shown maps are a combination of OpenStreetMaps and the image, which I generate like shown above. I use a tool called QGIS which allows me to combine different geo-referenced layers. QGIS allows me to manipulate these layers in color and transparancy etc.

pascalgross commented 4 years ago

Setting the unmapped points to nan is solving the whole problem. The rasterio library handles them correctly and the image is generated as wanted! Thank you!

JiaweiZhuang commented 4 years ago

Great, glad to hear that