These are chat archives for SHTOOLS/SHTOOLS

29th
Nov 2016
MMesch
@MMesch
Nov 29 2016 09:55
just mentioning, this is not a cylindrical projection if lat values are on the GLQ grid. We should interpolate to be on an even lat/lon grid. Then plot either without or without whatever library for projections ...
        """Plot the raw data using a simple cylindrical projection."""
        fig, ax = _plt.subplots(2, 1)
        ax.flat[0].imshow(self.data.real, origin='upper')
@MarkWieczorek , if the south-pole is undefined as in the DH grids, should we just interpolate (=average) the sourrounding points?
Mark Wieczorek
@MarkWieczorek
Nov 29 2016 10:03
Well, for DH grids it is cylindrical. If you really want to project the GLQ grids then I guess we would need to interpolate.
For the south pole, this has always been annoying. Driscoll and Healy use N by N grids, where the north pole is counted, but no the south pole. Whats worse, is that when you do the integrations in the spherical harmonic expansion, the north pole is downweighted to zero. For the plotting routines, averaging the last band would be good enough. But we might want to give an option in SHCoeffs.expand() to calculate the south pole value explicitly.
MMesch
@MMesch
Nov 29 2016 10:12
agree. I don't understand why they didn't choose the grid to be symmetric around the equator...
Mark Wieczorek
@MarkWieczorek
Nov 29 2016 10:13
I might have the fortran stop PR done tonight. If not, tomorrow morning.
MMesch
@MMesch
Nov 29 2016 10:17
good. I can work on a plot function today
so we should have everything ready very soon
Mark Wieczorek
@MarkWieczorek
Nov 29 2016 10:19
After this, I just want to implement expanding to a series of grid points, and that is it. That should be done early next week.
MMesch
@MMesch
Nov 29 2016 10:19
ok
we could even use this expand function for the plotting routines. It would be more accurate than using whatever other interpolator.
I don't know how fast/slow it is
Mark Wieczorek
@MarkWieczorek
Nov 29 2016 10:22
10s to 100s of times slower if you want to do all points. For plots, interpolation is fine with me, but I wouldn't be happy with giving the user the interpolated grid itself to work with unless this is all documented.
MMesch
@MMesch
Nov 29 2016 14:05
I added a plot function using cartopy now. Seems to work well...
The DH plot is weird.
MMesch
@MMesch
Nov 29 2016 14:30
I think there is a bug Mark:
In [47]: coeffs = pyshtools.SHCoeffs.from_zeros(10)

In [48]: coeffs.set_coeffs(1, 1, 0)

In [49]: grid = coeffs.expand(grid='DH')

In [50]: grid.lats()
Out[50]: 
array([ 90.        ,  81.81818182,  73.63636364,  65.45454545,
        57.27272727,  49.09090909,  40.90909091,  32.72727273,
        24.54545455,  16.36363636,   8.18181818,   0.        ,
        -8.18181818, -16.36363636, -24.54545455, -32.72727273,
       -40.90909091, -49.09090909, -57.27272727, -65.45454545,
       -73.63636364, -81.81818182])

In [51]: grid.data[:, 0]
Out[51]: 
array([ 1.73205081,  1.71442103,  1.66189058,  1.57552883,  1.45709386,
        1.30899666,  1.13425206,  0.93641736,  0.71951991,  0.4879751 ,
        0.24649653,  0.        , -0.24649653, -0.4879751 , -0.71951991,
       -0.93641736, -1.13425206, -1.30899666, -1.45709386, -1.57552883,
       -1.66189058, -1.71442103])

In [52]: grid.plot()
Out[52]: 
(<matplotlib.figure.Figure at 0x7fd1ab4495c0>,
 <matplotlib.axes._subplots.AxesSubplot at 0x7fd1ab449fd0>)
DHgrid_bug.png
sorry, it seems to be correct. l1m0 coefficient is symmetric
and you just have the additional point at the north pole
MMesch
@MMesch
Nov 29 2016 14:37
cartopy plotting works well now...
in case of the DH grid, we should explicitly calculate the South-Pole
MMesch
@MMesch
Nov 29 2016 14:42
also we could think about using interpolation to return high resolution grids of low order spectra, in case it is faster than a true expansion.
Mark Wieczorek
@MarkWieczorek
Nov 29 2016 17:22
I'm done with the Fortran stop PR. Please read my final comments, and then I will merge this tomorrow: SHTOOLS/SHTOOLS#59
Mark Wieczorek
@MarkWieczorek
Nov 29 2016 22:24
@MMesch what if instead of the keyword coordinates we used projection instead? I think that PlateCarre is cylindrical, and we could alias the two if necessary. Also, perhaps a projection=None would just plot the matrix (or perhaps default to cylindrical). Also, we should allow the method to pass generic keyword args to this somehow:
plt.axes(projection=ccrs.Orthographic(central_longitude=0.0,
 +                           central_latitude=80., globe=None))
Mark Wieczorek
@MarkWieczorek
Nov 29 2016 22:32
Also, for the normal imshow plots, maybe we should add interpolation=bicubic as the default, because the default of None looks horrible.
Also, for the interpolation you set
nlons_new, nlats_new = 400, 200
which is fine for the example you are using, but what would happen if the initial grid was larger than this?