Permalink
Browse files

longitude values are limited to the given range during the conversion

  • Loading branch information...
1 parent e332082 commit 1760142a5e5b77a380299468b1a59868a48436ce @leejjoon committed Feb 11, 2011
Showing with 164 additions and 7 deletions.
  1. +87 −4 lib/axes_wcs.py
  2. +77 −3 lib/wcs_helper.py
View
@@ -789,6 +789,10 @@ def set_ticklabel_type(self, labtyp1, labtyp2,
self.set_ticklabel2_type(labtyp2, **labtyp2_kwargs)
+ def set_lon_ref(self, ref):
+ #_fix_lon
+ self.projection.set_lon_ref(ref)
+
class GridHelperWcsSkyBase(GridHelperWcsBase):
_GRIDHELPER_CLASS=GridHelperCurveLinear
@@ -857,6 +861,83 @@ def get_display_coord_system(self):
return self._wcsgrid_display_coord_system
+ def _get_line_path_naive(self, axes, x, y):
+
+ from matplotlib.path import Path
+
+ if self.projection._lon_ref is None:
+ return Path(zip(x, y))
+
+ #tr = axes[0].transAux.inverted()
+ #rrdd = tr.transform(np.array([x, y]).transpose())
+
+ #rr, dd = rrdd[:,0], rrdd[:,1]
+ rr, dd = self.projection.toworld((x,y))
+ #print rr.astype("i"), rr1.astype("i")
+ rr += self.projection._lon_ref #self._lon_ref
+
+ rri = np.floor_divide(rr, 360.)
+ rrid = rri[1:] - rri[:-1]
+
+ codes = np.empty(rri.shape, dtype="i")
+ codes.fill(Path.LINETO)
+
+ i0 = 0
+ for i1, in np.transpose(rrid.nonzero()):
+ codes[i0] = Path.MOVETO
+ i0 = i1+1
+ codes[i0] = Path.MOVETO
+
+ return Path(zip(x, y), codes=codes)
+
+
+ def _get_line_path(self, axes, x, y):
+
+ from matplotlib.path import Path
+
+ if self.projection._lon_ref is None:
+ return Path(zip(x, y))
+
+ #rr, dd = self.projection.toworld((x,y))
+ #rrid = np.abs(rr[1:] - rr[:-1]) > 50.
+ x1, x2 = axes.get_xlim()
+ dx = abs(x2-x1)*.1
+ rrid = np.abs(x[1:] - x[:-1]) > dx
+ codes = np.empty(x.shape, dtype="i")
+ codes.fill(Path.LINETO)
+
+ i0 = 0
+ codes[i0] = Path.MOVETO
+ for i1, in np.transpose(rrid.nonzero()):
+ codes[i1+1] = Path.MOVETO
+
+ return Path(zip(x, y), codes=codes)
+
+
+ def new_floating_axis(self, nth_coord,
+ value,
+ axes=None,
+ axis_direction="bottom",
+ allsky=False
+ ):
+
+ #axisline = self._GRIDHELPER_CLASS.new_floating_axis( \
+ # self, nth_coord, value, axes, axis_direction)
+ axisline = super(GridHelperWcsSkyBase, self).new_floating_axis( \
+ nth_coord, value, axes, axis_direction)
+
+ if allsky:
+ h = axisline.get_helper()
+
+ # fixme
+ h._line_num_points = 1500
+ h._get_line_path = self._get_line_path
+
+ #axisline.line.set_clip_on(True)
+ #axisline.line.set_clip_box(axisline.axes.patch)
+ return axisline
+
+
class GridHelperWcsSimple(GridHelperWcsBase, GridHelperCurveLinear):
_GRIDHELPER_CLASS=GridHelperCurveLinear
@@ -949,7 +1030,7 @@ def __init__(self, wcs, orig_coord=None,
GridHelperWcs = GridHelperWcsSky
-class GridHelperWcsFloating(GridHelperWcsBase, floating_axes.GridHelperCurveLinear):
+class GridHelperWcsFloating(GridHelperWcsSkyBase, floating_axes.GridHelperCurveLinear):
def __init__(self, wcs, extremes,
orig_coord=None,
grid_locator1=None,
@@ -958,7 +1039,10 @@ def __init__(self, wcs, extremes,
tick_formatter2=None,
):
- GridHelperWcsBase.__init__(self, wcs, orig_coord)
+ #GridHelperWcsBase.__init__(self, wcs, orig_coord)
+
+ orig_coord, axis_nums = None, None
+ GridHelperWcsSkyBase.__init__(self, wcs, orig_coord, axis_nums)
if grid_locator1 is None:
@@ -1250,7 +1334,6 @@ def _init_parasites(self):
self._wcsgrid_wcsaxes = {0:ax}
self.parasites.append(ax)
-
def __getitem__(self, key):
# check if key is a valid coord_sys instance
@@ -1412,6 +1495,7 @@ def _decorate_default_label(self, label1, ticktyp1):
else:
if ticktyp1 == "absdeg":
label1 = label1 + r" [$^{\circ}$]"
+
return label1
def set_default_label(self, ticktyp1=None, ticktyp2=None):
@@ -1726,4 +1810,3 @@ def add_size_bar(self, length_pixel, label, loc, sep=5,
-
View
@@ -168,6 +168,45 @@ class ProjectionBase(object):
"""
A wrapper for kapteyn.projection or pywcs
"""
+
+ def __init__(self):
+ self._lon_ref = None
+
+ for i, ct in enumerate(self.ctypes):
+ if _pattern_ra.match(ct) or _pattern_glon.match(ct):
+ self._lon_axis = i
+ break
+ else:
+ self._lon_axis = None
+
+ def get_lon(self):
+ pass
+
+
+ def set_lon_ref(self, ref):
+ self._lon_ref = ref
+
+ def _fix_lon(self, lon, lon_ref=None):
+ """
+ transform lon into values in range [self._lon_ref, self._lon_ref+360]
+ """
+ if lon_ref is None:
+ lon_ref = self._lon_ref
+
+ if self._lon_ref is not None:
+ lon_ = lon - self._lon_ref
+ lon2 = lon_ - 360*np.floor_divide(lon_, 360.)
+ lon2[lon_==360] = 360
+ return lon2 + self._lon_ref
+ else:
+ return lon
+
+ def fix_lon(self, lon_lat):
+ """
+ """
+ if self._lon_axis is not None:
+ lon_lat[self._lon_axis] = self._fix_lon(lon_lat[self._lon_axis])
+
# def _get_ctypes(self):
# pass
@@ -195,6 +234,7 @@ class ProjectionKapteyn(ProjectionBase):
A wrapper for kapteyn.projection
"""
def __init__(self, header):
+ ProjectionBase.__init__(self)
if isinstance(header, pyfits.Header):
self._proj = kapteyn.wcs.Projection(header)
else:
@@ -250,6 +290,8 @@ def __init__(self, header):
else:
self._pywcs = header
+ ProjectionBase.__init__(self)
+
def _get_ctypes(self):
return tuple(self._pywcs.wcs.ctype)
@@ -268,13 +310,29 @@ def _get_naxis(self):
def topixel(self, xy):
""" 1, 1 base """
- xy2 = self._pywcs.wcs_sky2pix(np.asarray(xy).T, 1)
- return xy2.T
+
+ lon_lat = np.array(xy)
+
+ self.fix_lon(lon_lat)
+
+ xy1 = lon_lat.transpose()
+
+ # somehow, wcs_sky2pix does not work for some cases
+ xy21 = [self._pywcs.wcs_sky2pix([xy11], 1)[0] for xy11 in xy1]
+ #xy21 = self._pywcs.wcs_sky2pix(xy1, 1)
+
+ xy2 = np.array(xy21).transpose()
+ return xy2
def toworld(self, xy):
""" 1, 1 base """
xy2 = self._pywcs.wcs_pix2sky(np.asarray(xy).T, 1)
- return xy2.T
+
+ lon_lat = xy2.T
+ # fixme
+ self.fix_lon(lon_lat)
+
+ return lon_lat
#def sub(self, axes):
# wcs = self._pywcs.sub(axes=axes)
@@ -301,6 +359,7 @@ def __init__(self, proj, axis_nums_to_keep, ref_pixel):
_ref_world0 = np.asarray(proj.toworld(_ref_pixel0))
self._ref_world = _ref_world0.reshape((len(ref_pixel),))
+ ProjectionBase.__init__(self)
# for n in range(proj.naxis):
@@ -377,6 +436,9 @@ def toworld(self, xy):
#xyz2r = [d for (i, d) in enumerate(xyz2) if i in self._axis_nums_to_keep]
xyz2r = [xyz2[i] for i in self._axis_nums_to_keep]
+ # fixme
+ #xyz2r[0] = self._fix_lon(xyz2r[0])
+
return xyz2r
@@ -386,6 +448,7 @@ class ProjectionPywcs(ProjectionBase):
A wrapper for pywcs
"""
def __init__(self, header):
+ ProjectionBase.__init__(self)
if isinstance(header, pyfits.Header):
self._pywcs = pywcs.WCS(header=header)
else:
@@ -423,6 +486,7 @@ class ProjectionSimple(ProjectionBase):
A wrapper for pywcs
"""
def __init__(self, header):
+ ProjectionBase.__init__(self)
self._pywcs = pywcs.WCS(header=header)
self._simple_init(header)
@@ -454,6 +518,9 @@ def _simple_init(self, header):
def _simple_to_pixel(self, lon, lat):
lon, lat = np.asarray(lon), np.asarray(lat)
+
+ lon = self._fix_lon(lon)
+
x = (lon - self.crval1)/self.cdelt1*self.cos_phi + self.crpix1
y = (lat - self.crval2)/self.cdelt2 + self.crpix2
@@ -464,6 +531,8 @@ def _simple_to_world(self, x, y):
lon = (x - self.crpix1)*self.cdelt1/self.cos_phi + self.crval1
lat = (y - self.crpix2)*self.cdelt2 + self.crval2
+ lon = self._fix_lon(lon)
+
return lon, lat
def topixel(self, xy):
@@ -491,6 +560,11 @@ def sub(self, axis_nums):
if _pywcs_installed:
+# def ProjectionPywcsNdSimple(header):
+# if header["ctype1"].lower().endswith("car") and \
+# header["ctype2"].lower().endswith("car"):
+# return ProjectionSimple
+# return ProjectionPywcsNd
ProjectionDefault = ProjectionPywcsNd
else:
ProjectionDefault = ProjectionKapteyn

0 comments on commit 1760142

Please sign in to comment.