BUG #1886: Polar projection does not change pole.

Packages/vcs/vcs/template.py

@@ -1595,7 +1595,7 @@ def plot(self, x, slab, gm, bg=0, min=None,
             Y = slab.getAxis(-2)
 
         wc2 = vcs.utils.getworldcoordinates(gm, X, Y)
-        wc2 = kargs.get("dataset_bounds", wc2)
+        wc2 = kargs.get("plotting_dataset_bounds", wc2)
 
         # Do the boxes and lines
         for tp in ["box", "line"]:

Packages/vcs/vcs/vcs2vtk.py

@@ -251,11 +251,18 @@ def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None,
     if geo is None:
         bounds = pts.GetBounds()
         xm, xM, ym, yM = [bounds[0], bounds[1], bounds[2], bounds[3]]
-        # We don't use the zooming feature (gm.datawc) for geographic
-        # projections. We use wrapped coordinates for doing the projection
-        # such that parameters such that central meridian are set correctly
+        # We use zooming feature (gm.datawc) for linear and polar projections.
+        # We use wrapped coordinates for doing the projection
+        # such that parameters like the central meridian are set correctly.
+        if (gm.g_name == 'Gfm'):
+            # axes are not lon/lat for meshfill
+            wc = [gm.datawc_x1, gm.datawc_x2, gm.datawc_y1, gm.datawc_y2]
+        else:
+            wc = vcs.utils.getworldcoordinates(gm,
+                                               data1.getAxis(-1),
+                                               data1.getAxis(-2))
         geo, geopts = project(pts, projection, getWrappedBounds(
-            [gm.datawc_x1, gm.datawc_x2, gm.datawc_y1, gm.datawc_y2], [xm, xM, ym, yM], wrap))
+            wc, [xm, xM, ym, yM], wrap))
         pts = geopts
     # Sets the vertices into the grid
     if grid is None:
@@ -497,12 +504,18 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None):
     projection = vcs.elements["projection"][gm.projection]
     if grid is None:
         vg.SetPoints(pts)
-        # Even if we don't use the zooming feature (gm.datawc) for geographic
-        # projections
-        # we use plotting coordinates for doing the projection
+        # We use the zooming feature for linear and polar projections
+        # We use plotting coordinates for doing the projection
         # such that parameters such that central meridian are set correctly
+        if (gm.g_name == 'Gfm'):
+            # axes are not lon/lat for meshfill
+            wc = [gm.datawc_x1, gm.datawc_x2, gm.datawc_y1, gm.datawc_y2]
+        else:
+            wc = vcs.utils.getworldcoordinates(gm,
+                                               data1.getAxis(-1),
+                                               data1.getAxis(-2))
         geo, geopts = project(pts, projection, getWrappedBounds(
-            [gm.datawc_x1, gm.datawc_x2, gm.datawc_y1, gm.datawc_y2], [xm, xM, ym, yM], wrap))
+            wc, [xm, xM, ym, yM], wrap))
         # Sets the vertics into the grid
         vg.SetPoints(geopts)
     else:
@@ -603,10 +616,7 @@ def apply_proj_parameters(pd, projection, x1, x2, y1, y2):
     projName = pname
     pd.SetName(projName)
     if projection.type == "polar (non gctp)":
-        minY = min(y1, y2)
-        maxY = max(y1, y2)
-        if ((minY + 90.0) <= (90.0 - maxY)):
-            # minY is closer to -90 than maxY to 90
+        if (y1 < y2):
             pd.SetOptionalParameter("lat_0", "-90.")
             pd.SetCentralMeridian(x1)
         else:
@@ -1875,8 +1885,11 @@ def vtkIterate(iterator):
         obj = iterator.GetNextItem()
 
 
-# Return gmbounds if gmbounds are different than 1.e20
 def getPlottingBounds(gmbounds, databounds, geo):
+    """
+    Returns databounds for geographic projection
+    else returns gmbounds if gmbounds are different than 1.e20
+    """
     x1gm, x2gm, y1gm, y2gm = gmbounds[:4]
     x1, x2, y1, y2 = databounds[:4]
     if geo:
@@ -1889,10 +1902,12 @@ def getPlottingBounds(gmbounds, databounds, geo):
     return [x1, x2, y1, y2]
 
 
-# transforms [v1,v2] and returns it
-# such that it is in the same order
-# and has the same middle interval as [gm1, gm2]
 def switchAndTranslate(gm1, gm2, v1, v2, wrapModulo):
+    """
+    Transforms [v1,v2] and returns it
+    such that it is in the same order
+    and has the same middle interval as [gm1, gm2]
+    """
     assert(v1 < v2)
     # keep the same middle of the interval
     if (wrapModulo):
@@ -1907,15 +1922,12 @@ def switchAndTranslate(gm1, gm2, v1, v2, wrapModulo):
     return [v1, v2]
 
 
-# Returns bounds with the same interval size as databounds
-# but in the same order and with the same middle interval
-# as gmbounds. The middle and the order are used for
-# plotting. wrapModule has YWrap, XWrap in degrees, 0 means no wrap
 def getWrappedBounds(gmbounds, databounds, wrapModulo):
-    """ Returns the same interval as databounds but it
-    matches the order and also it keeps the same center interval as gmbounds
-    So for instance if databounds is -40, 320 and gmbounds is -180, 180
-    this function returns
+    """
+    Returns bounds with the same interval size as databounds
+    but in the same order and with the same middle interval
+    as gmbounds. The middle and the order are used for
+    plotting. wrapModule has YWrap, XWrap in degrees, 0 means no wrap
     """
     x1gm, x2gm, y1gm, y2gm = gmbounds[:4]
     x1, x2, y1, y2 = databounds[:4]

testing/vcs/CMakeLists.txt

@@ -11,6 +11,11 @@ cdat_add_test(test_vcs_bad_png_path
   "${PYTHON_EXECUTABLE}"
   ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_bad_png_path.py
 )
+cdat_add_test(test_vcs_boxfill_polar
+  "${PYTHON_EXECUTABLE}"
+  ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_polar.py
+  "${BASELINE_DIR}/test_vcs_boxfill_polar.png"
+  )
 cdat_add_test(test_vcs_create_get
   "${PYTHON_EXECUTABLE}"
   ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_create_get.py

testing/vcs/test_vcs_boxfill_polar.py

@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+import cdms2
+import os
+import sys
+import vcs
+
+pth = os.path.join(os.path.dirname(__file__), "..")
+sys.path.append(pth)
+import checkimage
+
+f = cdms2.open(vcs.sample_data + "/clt.nc")
+a=f("clt")
+
+x=vcs.init()
+x.setantialiasing(0)
+x.drawlogooff()
+x.setbgoutputdimensions(1200, 900, units="pixels")
+
+p=x.getprojection("polar")
+b=x.createboxfill()
+b.projection=p
+#b.datawc_y1 = 90
+#b.datawc_y2 = -90
+
+x.setbgoutputdimensions(1200,1091,units="pixels")
+x.plot(a(latitude=(90,-90)), b, bg=1)
+
+fileName = os.path.basename(__file__)
+fileName = os.path.splitext(fileName)[0]
+fileName += '.png'
+x.png(fileName)
+ret = checkimage.check_result_image(fileName, sys.argv[1], checkimage.defaultThreshold)
+sys.exit(ret)