Fix vector plots

CMake/cmake_modules/FindMatplotlib.cmake

@@ -4,7 +4,7 @@ if ( NOT ${res} EQUAL 0 )
  message("Could not detect matplotlib, ${PYTHON_EXECUTABLE}")
  set( Matplotlib_FOUND 0 )
 else ( NOT ${res} EQUAL 0 )
- execute_process(COMMAND "@PYTHON_EXECUTABLE@" -c "import matplotlib;print matplotlib.__version__" RESULT_VARIABLE res OUTPUT_VARIABLE out ERROR_VARIABLE err)
+ execute_process(COMMAND ${PYTHON_EXECUTABLE} -c "import matplotlib;print matplotlib.__version__" RESULT_VARIABLE res OUTPUT_VARIABLE out ERROR_VARIABLE err)
  if ( ${out} VERSION_LESS ${MATPLOTLIB_MAJOR_MIN}.${MATPLOTLIB_MINOR_MIN}.${MATPLOTLIB_PATCH_MIN} )
    message("You need matplotlib version ${MATPLOTLIB_MAJOR_MIN}.${MATPLOTLIB_MINOR_MIN}.${MATPLOTLIB_PATCH_MIN} at minimum")
    set( Matplotlib_FOUND 0 )

Packages/vcs/Lib/VTKPlots.py

@@ -1155,7 +1155,6 @@ def fitToViewport(self, Actor, vp, wc=None, geo=None, priority=None,
                         pt.InsertPoint(NGridCover, x, y, 0)
                         NGridCover += 1
                 pts = vtk.vtkPoints()
-                # pts.SetNumberOfPoints(Npts*Npts)
                 geo.TransformPoints(pt, pts)
                 b = pts.GetBounds()
                 xm, xM, ym, yM = b[:4]
@@ -1169,6 +1168,7 @@ def fitToViewport(self, Actor, vp, wc=None, geo=None, priority=None,
                     Yrg2[1] = max(Yrg2[1], yM)
                 Xrg = Xrg2
                 Yrg = Yrg2
+
             wRatio = float(sc[0]) / float(sc[1])
             dRatio = (Xrg[1] - Xrg[0]) / (Yrg[1] - Yrg[0])
             vRatio = float(vp[1] - vp[0]) / float(vp[3] - vp[2])

Packages/vcs/Lib/vcs2vtk.py

@@ -93,10 +93,6 @@ def putMaskOnVTKGrid(data, grid, actorColor=None, cellData=True, deep=True):
                 grid2.GetPointData().RemoveArray(
                     vtk.vtkDataSetAttributes.GhostArrayName())
                 grid2.GetPointData().SetScalars(imsk)
-                # grid2.SetCellVisibilityArray(imsk)
-                # p2c = vtk.vtkPointDataToCellData()
-                # p2c.SetInputData(grid2)
-                # geoFilter.SetInputConnection(p2c.GetOutputPort())
                 geoFilter.SetInputData(grid2)
                 lut.SetTableValue(0, r / 100., g / 100., b / 100., 1.)
                 lut.SetTableValue(1, r / 100., g / 100., b / 100., 1.)
@@ -137,9 +133,7 @@ def handleProjectionEdgeCases(projection, data):
     # For mercator projection, latitude values of -90 or 90
     # transformation result in infinity values. We chose -85, 85
     # as that's the typical limit used by the community.
-
     pname = projDict.get(projection._type, projection.type)
-
     if (pname.lower() == "merc"):
         lat = data.getLatitude()[:]
         # Reverse the latitudes incase the starting latitude is greater
@@ -150,13 +144,16 @@ def handleProjectionEdgeCases(projection, data):
     return data
 
 
-def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None):
+def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None,
+                    data2=None):
     continents = False
     projection = vcs.elements["projection"][gm.projection]
     xm, xM, ym, yM = None, None, None, None
     useStructuredGrid = True
 
     data1 = handleProjectionEdgeCases(projection, data1)
+    if data2 is not None:
+        data2 = handleProjectionEdgeCases(projection, data2)
 
     try:
         g = data1.getGrid()
@@ -216,14 +213,15 @@ def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None):
     xm, xM, ym, yM = getRange(gm, xm, xM, ym, yM)
     if geo is None:
         geo, geopts = project(pts, projection, [xm, xM, ym, yM])
+        pts = geopts
     # Sets the vertices into the grid
     if grid is None:
         if useStructuredGrid:
             vg = vtk.vtkStructuredGrid()
             vg.SetDimensions(data1.shape[1], data1.shape[0], 1)
         else:
             vg = vtk.vtkUnstructuredGrid()
-        vg.SetPoints(geopts)
+        vg.SetPoints(pts)
     else:
         vg = grid
     out = {"vtk_backend_grid": vg,
@@ -234,7 +232,8 @@ def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None):
            "continents": continents,
            "wrap": wrap,
            "geo": geo,
-           "data": data1
+           "data": data1,
+           "data2": data2
            }
     return out
 
@@ -546,6 +545,41 @@ def prepContinents(fnm):
     return poly
 
 
+def projectArray(w, projection, wc, geo=None):
+    xm, xM, ym, yM = wc
+    if isinstance(projection, (str, unicode)):
+        projection = vcs.elements["projection"][projection]
+    if projection.type == "linear":
+        return None, w
+
+    if geo is None:
+        geo = vtk.vtkGeoTransform()
+        ps = vtk.vtkGeoProjection()
+        pd = vtk.vtkGeoProjection()
+
+        pname = projDict.get(projection._type, projection.type)
+        projName = pname
+        pd.SetName(projName)
+
+        if projection.type == "polar (non gctp)":
+            if ym < yM:
+                pd.SetOptionalParameter("lat_0", "-90.")
+                pd.SetCentralMeridian(xm)
+            else:
+                pd.SetOptionalParameter("lat_0", "90.")
+                pd.SetCentralMeridian(xm + 180.)
+        else:
+            setProjectionParameters(pd, projection)
+        geo.SetSourceProjection(ps)
+        geo.SetDestinationProjection(pd)
+
+    for i in range(0, w.GetNumberOfTuples()):
+        tuple = [0, 0, 0]
+        w.GetTupleValue(i, tuple)
+        geo.TransformPoint(tuple, tuple)
+        w.SetTupleValue(i, tuple)
+
+
 # Geo projection
 def project(pts, projection, wc, geo=None):
     xm, xM, ym, yM = wc

Packages/vcs/Lib/vcsvtk/vectorpipeline.py

@@ -4,6 +4,8 @@
 from vcs import vcs2vtk
 import vtk
 
+import math
+
 
 class VectorPipeline(Pipeline):
 
@@ -16,6 +18,7 @@ def plot(self, data1, data2, tmpl, grid, transform):
         """Overrides baseclass implementation."""
         # Preserve time and z axis for plotting these inof in rendertemplate
         geo = None  # to make flake8 happy
+        projection = vcs.elements["projection"][self._gm.projection]
         returned = {}
         taxis = data1.getTime()
         if data1.ndim > 2:
@@ -27,14 +30,50 @@ def plot(self, data1, data2, tmpl, grid, transform):
         data1 = self._context().trimData2D(data1)
         data2 = self._context().trimData2D(data2)
 
+        scale = 1.0
+        lat = None
+        lon = None
+
+        latAccessor = data1.getLatitude()
+        lonAccesrsor = data1.getLongitude()
+        if latAccessor:
+            lat = latAccessor[:]
+        if lonAccesrsor:
+            lon = lonAccesrsor[:]
+
+        if None not in (projection, lat, lon):
+            scale = (lat.max() - lat.min()) * (lon.max() - lon.min())
+
         gridGenDict = vcs2vtk.genGridOnPoints(data1, self._gm, deep=False, grid=grid,
-                                              geo=transform)
+                                              geo=transform, data2=data2)
+
+        data1 = gridGenDict["data"]
+        data2 = gridGenDict["data2"]
+        geo = gridGenDict["geo"]
+
         for k in ['vtk_backend_grid', 'xm', 'xM', 'ym', 'yM', 'continents',
                   'wrap', 'geo']:
             exec("%s = gridGenDict['%s']" % (k, k))
         grid = gridGenDict['vtk_backend_grid']
         self._dataWrapModulo = gridGenDict['wrap']
 
+        if geo is not None:
+            newv = vtk.vtkDoubleArray()
+            newv.SetNumberOfComponents(3)
+            newv.InsertTupleValue(0, [lon.min(), lat.min(),  0])
+            newv.InsertTupleValue(1, [lon.max(), lat.max(),  0])
+
+            vcs2vtk.projectArray(newv, projection,
+                                 [gridGenDict['xm'], gridGenDict['xM'],
+                                  gridGenDict['ym'], gridGenDict['yM']])
+            dimMin = [0, 0, 0]
+            dimMax = [0, 0, 0]
+            newv.GetTupleValue(0, dimMin)
+            newv.GetTupleValue(1, dimMax)
+            scale = (dimMax[1] - dimMin[1]) * (dimMax[0] - dimMin[0])/scale
+        else:
+            scale = 1.0
+
         returned["vtk_backend_grid"] = grid
         returned["vtk_backend_geo"] = geo
         missingMapper = vcs2vtk.putMaskOnVTKGrid(data1, grid, None, False,
@@ -68,6 +107,7 @@ def plot(self, data1, data2, tmpl, grid, transform):
 
         arrow = vtk.vtkGlyphSource2D()
         arrow.SetGlyphTypeToArrow()
+        arrow.SetOutputPointsPrecision(vtk.vtkAlgorithm.DOUBLE_PRECISION)
         arrow.FilledOff()
 
         glyphFilter = vtk.vtkGlyph2D()
@@ -81,7 +121,7 @@ def plot(self, data1, data2, tmpl, grid, transform):
 
         # Scale to vector magnitude:
         glyphFilter.SetScaleModeToScaleByVector()
-        glyphFilter.SetScaleFactor(2. * self._gm.scale)
+        glyphFilter.SetScaleFactor(math.sqrt(scale) * 2.0 * self._gm.scale)
 
         # These are some unfortunately named methods. It does *not* clamp the
         # scale range to [min, max], but rather remaps the range
@@ -90,7 +130,11 @@ def plot(self, data1, data2, tmpl, grid, transform):
         glyphFilter.SetRange(0.01, 1.0)
 
         mapper = vtk.vtkPolyDataMapper()
-        mapper.SetInputConnection(glyphFilter.GetOutputPort())
+
+        glyphFilter.Update()
+        data = glyphFilter.GetOutput()
+
+        mapper.SetInputData(data)
         act = vtk.vtkActor()
         act.SetMapper(mapper)
 
@@ -100,21 +144,26 @@ def plot(self, data1, data2, tmpl, grid, transform):
 
         x1, x2, y1, y2 = vcs.utils.getworldcoordinates(self._gm, data1.getAxis(-1),
                                                        data1.getAxis(-2))
+        if geo is None:
+            wc = [x1, x2, y1, y2]
+        else:
+            wc = None
+
+        # TODO: doWrap is broken for vectors
+        # act = vcs2vtk.doWrap(act, [x1, x2, y1, y2], self._dataWrapModulo)
 
-        act = vcs2vtk.doWrap(act, [x1, x2, y1, y2], self._dataWrapModulo)
         self._context().fitToViewport(act, [tmpl.data.x1, tmpl.data.x2,
                                             tmpl.data.y1, tmpl.data.y2],
-                                      [x1, x2, y1, y2],
+                                      wc=wc,
                                       priority=tmpl.data.priority,
                                       create_renderer=True)
 
-        returned.update(
-            self._context().renderTemplate(tmpl, data1, self._gm, taxis, zaxis))
+        returned.update(self._context().renderTemplate(tmpl, data1,
+                        self._gm, taxis, zaxis))
 
         if self._context().canvas._continents is None:
             continents = False
         if continents:
-            projection = vcs.elements["projection"][self._gm.projection]
             self._context().plotContinents(x1, x2, y1, y2, projection,
                                            self._dataWrapModulo, tmpl)
 

testing/vcs/CMakeLists.txt

@@ -43,7 +43,7 @@ if (NOT EXISTS /etc/redhat-release)
     cdat_add_test(vcs_test_antialiasing
       "${PYTHON_EXECUTABLE}"
       ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_antialiasing.py
-    )
+)
 endif()
 cdat_add_test(vcs_verify_init
   "${PYTHON_EXECUTABLE}"
@@ -702,6 +702,11 @@ cdat_add_test(vcs_test_taylor_2_quads
         ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_patterns.py
         "${BASELINE_DIR}/test_vcs_patterns.png"
         )
+      cdat_add_test(vcs_test_vectors_robinson
+        "${PYTHON_EXECUTABLE}"
+        ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_vectors_robinson.py
+        "${BASELINE_DIR}/test_vcs_vectors_robinson.png"
+        )
     endif()
   endif()
 

testing/vcs/test_vcs_vectors_robinson.py

@@ -0,0 +1,24 @@
+import vcs, cdms2, numpy, os, sys
+src = sys.argv[1]
+pth = os.path.join(os.path.dirname(__file__),"..")
+sys.path.append(pth)
+import checkimage
+
+x = vcs.init()
+x = vcs.init()
+x.setantialiasing(0)
+x.drawlogooff()
+x.setbgoutputdimensions(1200, 1091, units="pixels")
+f = cdms2.open(os.path.join(vcs.sample_data, "clt.nc"))
+u = f("u")
+v = f("v")
+V = x.createvector()
+p = x.createprojection()
+p.type = "robinson"
+V.projection = p
+x.plot(u,v,V, bg=1)
+
+fnm= "test_vcs_vectors_robinson.png"
+x.png(fnm)
+ret = checkimage.check_result_image(fnm, src, checkimage.defaultThreshold)
+sys.exit(ret)