NCAR / GeoCAT

GeoCAT website
https://geocat.ucar.edu/
14 stars 12 forks source link

xlim issue when crossing the dateline #63

Closed pletzer closed 1 year ago

pletzer commented 1 year ago

Describe the bug

I'm struggling with plotting a field over a sub-region that straddles the data line (longitude = 180 deg). The problem, which seems to be related to calling ax.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]], proj), shows up as only the west side of the plot being displayed. The original field is defined for longitudes = 0, ... 360.

To Reproduce

Run the following code to reproduce the problem:

import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import geocat.viz as gv
import matplotlib.pylab as plt

def create_latlon_data(nlat, nlon):
    """
    Create an xarray field over the entire earth, longitudes = 0...360 deg
    :param nlat: number of latitude cells
    :param nlon: number of longitude cells
    :returns xarray field 
    """
    nlat1, nlon1 = nlat + 1, nlon + 1

    # create a data array with lat-lon cooordinates
    lat = np.linspace(-90., 90, nlat1)

    # create a lon coordinate, the last value (ie 360 deg) is missing
    dlon = 360/nlon
    lon = np.linspace(0., 360, nlon1)

    xx, yy = np.meshgrid(lon, lat)
    data = np.sin(xx*np.pi/180.) * np.cos(yy*np.pi/180.)

    da = xr.DataArray(data, coords=[lat, lon], \
                      dims=['latitude', 'longitude'], name='foo')

    return da

# create data 
da = create_latlon_data(180, 360)

# quick plot to show what the data look like
#da.plot.contourf()

# now we want to plot the field between longitude 150 -> 200. Try different 
# ways of specifying xlim
ylim = (-50, -20)
for xlim in (150, 200), (150, 200 - 360):

    fig = plt.figure(figsize=(10, 8))
    proj = ccrs.PlateCarree()
    ax = plt.axes(projection=proj)
    ax.coastlines(linewidths=0.5)
    ax.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]], proj)
    gv.set_axes_limits_and_ticks(ax, xlim=xlim, ylim=ylim,
                                 xticks=np.linspace(xlim[0], xlim[1], 4+1),
                                 yticks=np.linspace(ylim[0], ylim[1], 4+1))
    gv.add_lat_lon_ticklabels(ax)
    da.plot.contourf(ax=ax, transform=proj, vmin=-1, vmax=1, levels=11,
                           cbar_kwargs={
                               "orientation": "horizontal",
                               "shrink": 0.8
                               }
                )
    plt.title(f'xlim = {xlim}')
    plt.show()

Expected behavior The above code produces two plots using different xlim values. In both cases the displayed domain is wrong. For the first plot, xlim = (150, 200) no data are shown from 180W to 160W. For the second plot, the field is displayed from 160W to 150E. See below.

Screenshots

image image

Environment (please complete the following information):

Additional context

kafitzgerald commented 1 year ago

I'm sorry we didn't see this sooner.

It looks like you might need this helper function to deal with the cyclic coordinates: https://geocat-viz.readthedocs.io/en/latest/user_api/generated/geocat.viz.util.xr_add_cyclic_longitudes.html

There's also an example here that might help: https://geocat-examples.readthedocs.io/en/latest/gallery/Contours/NCL_conLev_4.html#sphx-glr-gallery-contours-ncl-conlev-4-py

Let us know if you're still having trouble or I misunderstood what you're looking to do.

kafitzgerald commented 1 year ago

I'm going to close this for now, but if you're still having trouble with this feel free to reopen.