pangeo-data / xESMF

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

Use regrided data using xESMF to be plot with Matplotlib and Cartopy #262

Open nuvolet opened 1 year ago

nuvolet commented 1 year ago

Greetings,

I am trying to figure out why I cannot have a complete plot of from a data array I regrided to another resolution.

The original data array looks like that:

>>> Orig_data
<xarray.DataArray 'VISdn' (lat: 192, lon: 288)>
dask.array<mean_agg-aggregate, shape=(192, 288), dtype=float32, chunksize=(192, 288), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 87.17 88.12 89.06 90.0
  * lon      (lon) float64 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8

While after using the libraries xESMF:

ds_out = xe.util.grid_2d(0.0, 360.0, 1.0, -90.0, 90.0, 1.0)
regrid_Orig_data = xe.Regridder(Orig_data, ds_out, 'bilinear', periodic=True)
Orig_new_reg = regrid_Orig_data(Orig_data)

I got that:

>>> Orig_new_reg 
<xarray.DataArray (y: 180, x: 360)>
dask.array<_regrid, shape=(180, 360), dtype=float32, chunksize=(180, 360), chunktype=numpy.ndarray>
Coordinates:
    lat      (y, x) float64 -89.5 -89.5 -89.5 -89.5 ... 89.5 89.5 89.5 89.5
    lon      (y, x) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
Dimensions without coordinates: y, x
Attributes:
    regrid_method:  bilinear

Then, when I plot the new data (Orig_new_reg ) using Matplotlib and Cartopy only shows my half of the data. Is this due to the new data says it is "Dimensions without coordinates: y, x"?. How I can solve that?

Thanks in advance

aulemahal commented 1 year ago

Is it only showing the half between 0 and 180 of longitude ? Could you maybe show us the plotting code with an image as example ?

I'm making the hypothesis that ESMF works in the -180-180 longitude system, not the 0-360 as you are using. Suggestion: replace your first line with:

ds_out = xe.util.grid_global(1, 1)

It does the same thing, is simpler and uses the -180-180 system.

huard commented 1 year ago

@nuvolet Is this resolved ?

nuvolet commented 1 year ago

Hi @aulemahal and @huard

I still cannot have a complete plot with data re-griided using xesmf. For instance the code I am using to plot the data I re-gridded looks like this:

ax1 = plt.subplot(221, projection=ccrs.Robinson()) ax1.set_global(); ax1.coastlines(); #ax1.gridlines(); ax1 = nudge_Aqua.plot(transform=ccrs.PlateCarree(), cmap='coolwarm', center=False, vmin=-0.5, vmax=0.5, levels=10, cbar_kwargs={'label': 'AOD[ ]', 'shrink': 0.85, 'format': "%2.1f"}) plt.title('CAM6 Nudge - MATCH Aqua', fontsize=8, loc='left')

Where nudge_Aqua is the difference from 2 datasets I previously re-gridded using xe.Regridder.

The plot I get looks like this:

Screenshot 2023-07-05 at 9 43 58 PM

huard commented 1 year ago

The only thing I can think of is that your input data has nans. Also, it might be easier to debug with a simple plt.imshow to rule out projection issues.

Luiz-Octavioo commented 11 months ago

Did someone manage to solve this problem? I have a similar issue. Here is the code. I would greatly appreciate some assistance. Thanks in advance

import xarray as xr
import xesmf as xe
import glob
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# Lat lon formatter
from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter
# Add cyclic point
from cartopy.util import add_cyclic_point

# Get files to open
path = 'Data_CMIP/tx90p*ssp245*'
files = glob.glob(path)

# # Open files
datasets = []
for f in files:
    ds = xr.open_dataset(f).drop('height')
    try:
        ds['time'] = ds.indexes['time'].to_datetimeindex()
    except AttributeError:
        pass
    datasets.append(ds)

# Create grid to interpolation
ds_out = xe.util.grid_2d(-180.0, 180.0, 1.0, -90.0, 90.0, 1.0)
for i in range(len(datasets)):
    regridder = xe.Regridder(datasets[i], ds_out, 'bilinear', periodic=True)
    datasets[i] = regridder(datasets[i], keep_attrs=True)

# Concatenate all datasets into one
ds = xr.concat(datasets, dim='model', join='override')

# Calculate the multi-model mean
ensemble_mean = ds.tx90pETCCDI.mean(dim='model')

# Plot the multi-model mean
fig, ax = plt.subplots(figsize=(10, 5), subplot_kw={'projection': ccrs.Robinson(central_longitude=-180)})
ax.coastlines()
ax.set_global()
# ax.add_feature(cfeature.BORDERS)

# Add cyclic point
data_cyclic, lon_cyclic = add_cyclic_point(ensemble_mean, coord=ensemble_mean.x.values)
# Contourf
cs = ax.contourf(lon_cyclic, ensemble_mean.y.values, data_cyclic[-1], 
                 transform=ccrs.PlateCarree(),
                 levels=40,
                cmap=cmap_temp_ncl, extend='both', robust=True)
# Colorbar
cbar = plt.colorbar(cs, orientation='horizontal', shrink=0.7, pad=0.05, extend='both', aspect=35)
huard commented 11 months ago

Hi Luiz,

Could you post a link to a sample of the data so I can try to reproduce the problem on my end ?