Map with polar stereographic has empty gap

[lee1043]

Map with polar stereo-graphic shows empty gap for fixed column which not caused by lack of data. Below plot is coming from one dataset and gap located at different longitude.

Is this bug? or is there any alternative way to fill this gap?

I have used following script:

import vcs
import cdutil

def main():
  f = cdms.open('/work/lee1043/cdat/sample_input/clt.nc')
  d = f('clt')
  plot(d, 1)

  region_selector = cdutil.region.domain(longitude=(-180.,180.))
  d2 = d(region_selector)
  plot(d2, 2)

def plot(d, n):
  x=vcs.init()
  iso = x.createisofill()
  p = vcs.createprojection()
  p.type = int('-3')
  iso.projection = p
  xtra = {}
  xtra['latitude'] = (90.0,0.0)
  d = d(**xtra)
  x.plot(d,iso)
  x.png('test'+str(n)+'.png')

main()```

[lee1043]

@doutriaux1 Have any idea or suggestion?

[doutriaux1]

@lee1043 you need to read extra data to close the gap longitude=(-180,182) should do.

[lee1043]

@doutriaux1 It sounds like a nice trick and helps to narrow the gap. But I will also need to use my region_selector for calculation, not only for the plot. I've tried various bound options such as ccb, ccn, ooe, oon, etc., but they were not so helpful.

[doutriaux1]

@lee1043 try to set datawc_X1/x2 to 0/360 (-180,180) it might help as well.

[doutriaux1]

@danlipsa any hints for this?

[lee1043]

@doutriaux1 Thanks but I have checked 0/360 or -180/180 is not helpful for this case.

[danlipsa]

@lee1043 @doutriaux1 Try (-180, 185). That will close the gap. The original data goes from (-180, 175) so a gap is appropriate. Not sure why (-180, 180) does not close it - I thought it should have. We'll have to look closer to see why.

[lee1043]

@danlipsa @doutriaux1 Is there any option that let the program know the data is cyclic so it can close the gap when plotting?

[doutriaux1]

hum... did youtrytweaking datwc_1/x2?

[lee1043]

@doutriaux1 I have tried adding datawc_x1 and datawc_x2 but it was not so helpful, still the gap is existing.

[lee1043]

Revisiting the issue with willingness of following up...

[doutriaux1]

@lee1043 does extending the data at plot time work?

x = vcs.init()
x.plot(data(longitude=(-180,185)))

[lee1043]

@doutriaux1 I know that could be a temporary solution but think it would be nice if that could be more generalized without having special treatment.

[lee1043]

@doutriaux1 this is interesting... I learned the gap only appears for isofill, but not for boxfill.

• boxfill:
  

• isofill:
  

• Plotting code: plot_test_simple.py.txt

import cdms2, vcs
f=cdms2.open('NAM_psl_EOF1_DJF_obs_1900-2005.nc')
d=f('eof')
d.getLongitude()

x=vcs.init()
gm=vcs.createboxfill()
gm.projection='polar'
x.plot(d,gm)
x.png('test_simple1.png')

x.clear()
gm=vcs.createisofill()
gm.projection='polar'
x.plot(d,gm)
x.png('test_simple2.png')

• Input data: NAM_psl_EOF1_DJF_obs_1900-2005.nc.txt
  (rename to .nc, removing .txt at the end)

[lee1043]

Another question: Can the stitches along the 0 lon (probably latitude tick label?) be removed?

[doutriaux1]

@lee1043 yes turn them off via template. template.ytic1.priority=0 or something like that.

[lee1043]

@doutriaux1 thanks for the ytic1 hint. I learned both 'ytic1' and 'ytic2' should be considered as priority=0 to completely remove the stitches. What about the isofill thing? Do you have any idea? Thanks!

[doutriaux1]

I need to look into that one, I think I need to change isofill.

[jypeter]

It would indeed be nice to fix this gap problem automatically, rather than having to generate extra longitude values (even if cdms2 will nicely do this for you).

I've also had to generate extra longitudes for some non polar plots in order to avoid a white vertical band on the left or right side (don't remember which side) from time to time but I don't remember for which kind of data. Oh, you can actually see this in the Gallery:

• Problem with https://uvcdat.llnl.gov/examples/Clim_DJF_200mb_wnd_ts_gph.html
• NO side problem with https://uvcdat.llnl.gov/examples/Clim_JJA_precip.html

One last note: in the examples above, you plot either the North or South pole by plotting data(latitude=(90, -90)) or data(latitude=(-90, 90)), but I'm not sure this trick is explained somewhere. At least I have not found an obvious polar tutorial for this on the web site

[lee1043]

Agreed to @jypeter that we need tutorial for the polar stereographic map projection.

Just for future reference, latitude can be set as data(latitude=(90,0)) or data(latitude=(-90,0)) for showing only northern or southern hemisphere.

[jypeter]

Good point @lee1043 ! We just need half of the data for a polar plot :-)

[doutriaux1]

@jypeter actually the "polar" plot can display all 180 degrees of latitude. I'm not sure why it was historically called a "polar" plot, in reality it is an "Albers Equidistant" projection (proj4: aeqd)

[lee1043]

@doutriaux1 Albers Equidistant is slightly different then polar projection. Difference is at where to set edge point of flat projection of globe. Albers uses virtual point vertically above polar pole while polar projection uses actual pole at surface as one side of projected flat plane. Sorry about the poor explaining in text, I may able to give better description in person when I back ;)

[doutriaux1]

@lee1043 agreed, I just think vcs polar projection has been miscalled since the beginning.