LSSTDESC / skyproj

Sky mapping with matplotlib and PROJ
BSD 3-Clause "New" or "Revised" License
10 stars 5 forks source link

Lambert azimuthal equal area projection centered at lat_0 = -90 #22

Open kadrlica opened 2 years ago

kadrlica commented 2 years ago

Given that we are designing with LSST in mind, I think that it would probably make sense for the default of LaeaSkyproj to be centered at lat_0 = -90.

It seems that skyproj has trouble rendering this

smap = skyproj.LaeaSkyproj(lat_0=-90)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [30], in <module>
----> 1 smap = skyproj.LaeaSkyproj(lat_0=-90)

File ~/software/skyproj/skyproj/skyproj.py:1296, in LaeaSkyproj.__init__(self, **kwargs)
   1295 def __init__(self, **kwargs):
-> 1296     super().__init__(projection_name='laea', **kwargs)

File ~/software/skyproj/skyproj/skyproj.py:125, in Skyproj.__init__(self, ax, projection_name, lon_0, gridlines, celestial, extent, longitude_ticks, autorescale, **kwargs)
    122 self._top_line = None
    123 self._bottom_line = None
--> 125 self._initialize_axes(extent)
    127 # Set up callbacks on axis zoom.
    128 self._xlc = self._ax.callbacks.connect('xlim_changed', self._change_axis)

File ~/software/skyproj/skyproj/skyproj.py:205, in Skyproj._initialize_axes(self, extent)
    202     self._aa.remove()
    203     self._aa = None
--> 205 self._set_axes_limits(extent, invert=False)
    206 self._create_axes(extent)
    207 self._set_axes_limits(extent, invert=self.do_celestial)

File ~/software/skyproj/skyproj/skyproj.py:329, in Skyproj._set_axes_limits(self, extent, invert)
    326 if len(extent) != 4:
    327     raise ValueError("Must specify extent as a 4-element array.")
--> 329 self._ax.set_extent(extent, lonlat=True)
    331 if self._aa is not None:
    332     self._aa.set_xlim(self._ax.get_xlim())

File ~/software/skyproj/skyproj/skyaxes.py:89, in SkyAxes.set_extent(self, extent, lonlat)
     86         y0 = np.min(xy[:, 1])
     87         y1 = np.max(xy[:, 1])
---> 89 self.set_xlim([x0, x1])
     90 self.set_ylim([y0, y1])

File ~/.conda/envs/skyproj/lib/python3.8/site-packages/matplotlib/axes/_base.py:3688, in _AxesBase.set_xlim(self, left, right, emit, auto, xmin, xmax)
   3685     right = xmax
   3687 self._process_unit_info([("x", (left, right))], convert=False)
-> 3688 left = self._validate_converted_limits(left, self.convert_xunits)
   3689 right = self._validate_converted_limits(right, self.convert_xunits)
   3691 if left is None or right is None:
   3692     # Axes init calls set_xlim(0, 1) before get_xlim() can be called,
   3693     # so only grab the limits if we really need them.

File ~/.conda/envs/skyproj/lib/python3.8/site-packages/matplotlib/axes/_base.py:3605, in _AxesBase._validate_converted_limits(self, limit, convert)
   3602 converted_limit = convert(limit)
   3603 if (isinstance(converted_limit, Real)
   3604         and not np.isfinite(converted_limit)):
-> 3605     raise ValueError("Axis limits cannot be NaN or Inf")
   3606 return converted_limit

ValueError: Axis limits cannot be NaN or Inf
posx and posy should be finite values
kadrlica commented 2 years ago

Oh, and for default lat_0=0 maps, it looks like full_sky isn't working and the map boundary is being drawn incorrectly.

image