Last active
December 21, 2020 18:06
-
-
Save matt-long/50433da346da8ac17cde926eec90a87c to your computer and use it in GitHub Desktop.
Plotting the POP model
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
import xarray as xr | |
def pop_add_cyclic(ds): | |
nj = ds.TLAT.shape[0] | |
ni = ds.TLONG.shape[1] | |
xL = int(ni/2 - 1) | |
xR = int(xL + ni) | |
tlon = ds.TLONG.data | |
tlat = ds.TLAT.data | |
tlon = np.where(np.greater_equal(tlon, min(tlon[:,0])), tlon-360., tlon) | |
lon = np.concatenate((tlon, tlon + 360.), 1) | |
lon = lon[:, xL:xR] | |
if ni == 320: | |
lon[367:-3, 0] = lon[367:-3, 0] + 360. | |
lon = lon - 360. | |
lon = np.hstack((lon, lon[:, 0:1] + 360.)) | |
if ni == 320: | |
lon[367:, -1] = lon[367:, -1] - 360. | |
#-- trick cartopy into doing the right thing: | |
# it gets confused when the cyclic coords are identical | |
lon[:, 0] = lon[:, 0] - 1e-8 | |
#-- periodicity | |
lat = np.concatenate((tlat, tlat), 1) | |
lat = lat[:, xL:xR] | |
lat = np.hstack((lat, lat[:,0:1])) | |
TLAT = xr.DataArray(lat, dims=('nlat', 'nlon')) | |
TLONG = xr.DataArray(lon, dims=('nlat', 'nlon')) | |
dso = xr.Dataset({'TLAT': TLAT, 'TLONG': TLONG}) | |
# copy vars | |
varlist = [v for v in ds.data_vars if v not in ['TLAT', 'TLONG']] | |
for v in varlist: | |
v_dims = ds[v].dims | |
if not ('nlat' in v_dims and 'nlon' in v_dims): | |
dso[v] = ds[v] | |
else: | |
# determine and sort other dimensions | |
other_dims = set(v_dims) - {'nlat', 'nlon'} | |
other_dims = tuple([d for d in v_dims if d in other_dims]) | |
lon_dim = ds[v].dims.index('nlon') | |
field = ds[v].data | |
field = np.concatenate((field, field), lon_dim) | |
field = field[..., :, xL:xR] | |
field = np.concatenate((field, field[..., :, 0:1]), lon_dim) | |
dso[v] = xr.DataArray(field, dims=other_dims+('nlat', 'nlon'), | |
attrs=ds[v].attrs) | |
# copy coords | |
for v, da in ds.coords.items(): | |
if not ('nlat' in da.dims and 'nlon' in da.dims): | |
dso = dso.assign_coords(**{v: da}) | |
return dso |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hi @matt-long:
I am using your util.py to plot sea surface temperature on POP gx1v6 grids. There are strange horizontal lines at the north pole. Please see the attached figure. My jupyter-notebook is available at Could you please let me know how to get rid of the horizontal lines near the north pole. Thanks so much for your help.
Regards!
Dapeng