Fixing vector plots not working for mercator and robinson projection

Packages/vcs/Lib/VTKPlots.py

@@ -640,6 +640,10 @@ def plotContinents(self, x1, x2, y1, y2, projection, wrap, tmpl):
         else:
             geo = None
 
+        writer = vtk.vtkXMLPolyDataWriter()
+        writer.SetFileName("poly.vtk")
+        writer.SetInputData(contData)
+        writer.Write()
         self.fitToViewport(contActor,
                            [tmpl.data.x1, tmpl.data.x2,
                             tmpl.data.y1, tmpl.data.y2],
@@ -1132,6 +1136,7 @@ def fitToViewport(self, Actor, vp, wc=None, geo=None, priority=None,
             else:
                 flipX = False
 
+
             if geo is not None:
                 pt = vtk.vtkPoints()
                 Xrg2 = [1.e20, -1.e20]
@@ -1150,7 +1155,7 @@ 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)
+                #pts.SetNumberOfPoints(Npts*Npts)
                 geo.TransformPoints(pt, pts)
                 b = pts.GetBounds()
                 xm, xM, ym, yM = b[:4]
@@ -1164,6 +1169,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

@@ -53,7 +53,6 @@
 for i in range(len(projNames)):
     projDict[i] = projNames[i]
 
-
 def applyAttributesFromVCStmpl(tmpl, tmplattribute, txtobj=None):
     tatt = getattr(tmpl, tmplattribute)
     if txtobj is None:
@@ -93,10 +92,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 +132,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
@@ -149,14 +142,16 @@ def handleProjectionEdgeCases(projection, data):
         data = data(latitude=(max(-85, lat.min()), min(85, lat.max())))
     return data
 
-
-def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None):
+def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None,
+                    skipReprojection=False, 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()
@@ -207,23 +202,25 @@ def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None):
         m3 = numpy.concatenate((m3, z), axis=1)
         deep = True
         pts = vtk.vtkPoints()
+        pts.SetDataTypeToDouble()
         # Convert nupmy array to vtk ones
         ppV = numpy_to_vtk_wrapper(m3, deep=deep)
         pts.SetData(ppV)
     else:
         xm, xM, ym, yM, tmp, tmp2 = grid.GetPoints().GetBounds()
         vg = grid
     xm, xM, ym, yM = getRange(gm, xm, xM, ym, yM)
-    if geo is None:
+    if geo is None and not skipReprojection:
         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 +231,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
 
@@ -545,6 +543,38 @@ def prepContinents(fnm):
     vcsContinents[fnm] = poly
     return poly
 
+def projectArray(w, projection, geo=None):
+    if isinstance(projection, (str, unicode)):
+        projection = vcs.elements["projection"][projection]
+    if projection.type == "linear":
+        return None, pts
+
+    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):
@@ -573,6 +603,7 @@ def project(pts, projection, wc, geo=None):
         geo.SetSourceProjection(ps)
         geo.SetDestinationProjection(pd)
     geopts = vtk.vtkPoints()
+    geopts.SetDataTypeToDouble()
     geo.TransformPoints(pts, geopts)
     return geo, geopts
 

Packages/vcs/Lib/vcsvtk/vectorpipeline.py

@@ -4,6 +4,7 @@
 from vcs import vcs2vtk
 import vtk
 
+import math
 
 class VectorPipeline(Pipeline):
 
@@ -16,6 +17,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 +29,44 @@ def plot(self, data1, data2, tmpl, grid, transform):
         data1 = self._context().trimData2D(data1)
         data2 = self._context().trimData2D(data2)
 
+        scale = 1.0
+
+        lat = data1.getLatitude()[:]
+        lon = data1.getLongitude()[:]
+
+        if projection is not None:
+            scale = (lat.max() - lat.min()) *  (lon.max() - lon.min())
+
         gridGenDict = vcs2vtk.genGridOnPoints(data1, self._gm, deep=False, grid=grid,
-                                              geo=transform)
+                                              geo=transform, skipReprojection=False,
+                                              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=projection)
+            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 +100,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 +114,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 +123,12 @@ def plot(self, data1, data2, tmpl, grid, transform):
         glyphFilter.SetRange(0.01, 1.0)
 
         mapper = vtk.vtkPolyDataMapper()
-        mapper.SetInputConnection(glyphFilter.GetOutputPort())
+        writer = vtk.vtkXMLPolyDataWriter()
+
+        glyphFilter.Update()
+        data = glyphFilter.GetOutput()
+
+        mapper.SetInputData(data)
         act = vtk.vtkActor()
         act.SetMapper(mapper)
 
@@ -100,21 +138,20 @@ def plot(self, data1, data2, tmpl, grid, transform):
 
         x1, x2, y1, y2 = vcs.utils.getworldcoordinates(self._gm, data1.getAxis(-1),
                                                        data1.getAxis(-2))
-
-        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=None,
                                       priority=tmpl.data.priority,
                                       create_renderer=True)
 
+
         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)