Arctic projection plotting routine

[NoraLoose]

Does anybody have a python plotting routine that produces nice Arctic projections from CCSM POP2 model output, e.g., with cartopy's NorthPolarStereo or Orthographic projections?

My first attempt was:

import cartopy
import cartopy.crs as ccrs

fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
ax.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND)
ax.gridlines()

cs = ax.pcolormesh(ds.TLAT,ds.TLONG,ds.TEMP.isel(z_t=0,time=-1),transform=ccrs.PlateCarree())

but it seems that there are interpolation issues (probably related to the seam of the tripolar grid in the Arctic):

This is of course not an issue with pop-tools, but I am hoping that this is a good platform to ask if anyone has experience with this. 😃 Thanks!

[matt-long]

Check out this function, pop_add_cyclic:
https://gist.github.com/matt-long/50433da346da8ac17cde926eec90a87c

You may also need to mask (subset) data outside the domain.

[Eddebbar]

Hi @NoraLoose, try switching ds.TLAT & ds.TLONG? i.e.:
cs = ax.pcolormesh(ds.TLONG,ds.TLAT,ds.TEMP.isel(z_t=0,time=-1),transform=ccrs.PlateCarree())

and for mask, this usually works for me:
land = ax.add_feature(cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',linewidth=0.1, edgecolor='black', facecolor='grey'))

[NoraLoose]

Thanks so much for your quick responses!

@Eddebbar, getting ds.TLAT and ds.TLONG in the right order was definitely a good idea 😄 ! But it didn't fix my problem completely, and I am still seeing weird interpolation patterns and stripes (despite masking out the data outside of the domain, as recommended by @matt-long).

@matt-long, I tried the add_cyclic_point function from #48 on my POP model output. But it seems that it did something else than shown in your notebook:

dsp = add_cyclic_point(ds)

fig,axs = plt.subplots(1,2,figsize=(18,5))
ds.TLONG.where(ds.KMT > 0).plot(ax=axs[0])
dsp.TLONG.where(dsp.KMT > 0).plot(ax=axs[1])

Should the add_cyclic_point function work on my POP_tx0.1 grid, too?

[Eddebbar]

@NoraLoose, hmmm, not sure what that's from... do the interpolation patterns/stripes look something like this? (seems like it's contouring over NaN/land)

fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
ax.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND)
ax.gridlines()
cs = ax.pcolormesh(ds.TLONG,ds.TLAT,ds.TEMP.isel(z_t=0,time=-1),transform=ccrs.PlateCarree())

which seems to go away on my end with cartopy.feature.NaturalEarthFeature, so could be a grid specific issue...

fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
ax.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND)
ax.gridlines()
cs = ax.pcolormesh(ds.TLONG,ds.TLAT,ds.TEMP.isel(z_t=0,time=-1),transform=ccrs.PlateCarree())
land = ax.add_feature(cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',linewidth=0.1, edgecolor='black', facecolor='white'))

[NoraLoose]

Sorry for not following up on your responses - I had dropped the ball on this issue for a while.

@Eddebbar, thanks for your suggestion! Weirdly, this didn't solve the issue on my end. I ended up decreasing the latitudinal extent from [50, 90] to [82, 90] because that's all I needed for our tutorial. Not a general solution, but oh well.

Thanks again for your help, @Eddebbar and @matt-long! Please feel free to close this issue.