Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
nuvolet opened this issue May 8, 2023 · 6 comments
Open

Comments

@nuvolet
Copy link

nuvolet commented May 8, 2023

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
Copy link
Collaborator

aulemahal commented May 8, 2023

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
Copy link
Contributor

huard commented Jun 15, 2023

@nuvolet Is this resolved ?

@nuvolet
Copy link
Author

nuvolet commented Jul 6, 2023

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
Copy link
Contributor

huard commented Jul 6, 2023

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
Copy link

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
Copy link
Contributor

huard commented Jul 17, 2023

Hi Luiz,

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants