painfully trying to plot isofill with gengrid

vcs/Canvas.py

@@ -2973,7 +2973,7 @@ def __plot(self, arglist, keyargs):
         if inGrid is not None and arglist[0] is not None and\
                 isinstance(arglist[0], cdms2.avariable.AbstractVariable) and\
                 not isinstance(arglist[0].getGrid(), cdms2.grid.AbstractRectGrid) and\
-                arglist[3] not in ["meshfill", ]:
+                arglist[3] not in ["meshfill", "isofill" ]:
             raise RuntimeError("You are attempting to plot unstructured grid" +
                                "with a method that is not meshfill")
         # preprocessing for extra keyword (at-plotting-time options)
@@ -3926,7 +3926,7 @@ def set_convert_labels(copy_mthd, test=0):
                     arglist[1].setAxis(i, axes_changed2[i])
             # Check to make sure that you have at least 2 dimensions for the follow graphics methods
             # Flipping the order to avoid the tv not exist problem
-            if (arglist[3] in ['boxfill', 'isofill', 'isoline', 'vector']) and (
+            if (arglist[3] in ['boxfill', 'isoline', 'vector']) and (
                     len(arglist[0].shape) < 2):
                 raise vcsError(
                     'Invalid number of dimensions for %s' %

vcs/vcs2vtk.py

@@ -340,21 +340,27 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False,
             continents = True
             wrap = [0., 360.]
             if grid is None:
-                m = g.getMesh()
-                xm = m[:, 1].min()
-                xM = m[:, 1].max()
-                ym = m[:, 0].min()
-                yM = m[:, 0].max()
-                numberOfCells = m.shape[0]
-                # For vtk we need to reorder things
-                m2 = numpy.ascontiguousarray(numpy.transpose(m, (0, 2, 1)))
-                m2.resize((m2.shape[0] * m2.shape[1], m2.shape[2]))
-                m2 = m2[..., ::-1]
-                # here we add dummy levels, might want to reconsider converting
-                # "trimData" to "reOrderData" and use actual levels?
-                m3 = numpy.concatenate(
-                    (m2, numpy.zeros(
-                        (m2.shape[0], 1))), axis=1)
+                try:
+                    m = g.getMesh()
+                    xm = m[:, 1].min()
+                    xM = m[:, 1].max()
+                    ym = m[:, 0].min()
+                    yM = m[:, 0].max()
+                    numberOfCells = m.shape[0]
+                    # For vtk we need to reorder things
+                    m2 = numpy.ascontiguousarray(numpy.transpose(m, (0, 2, 1)))
+                    m2.resize((m2.shape[0] * m2.shape[1], m2.shape[2]))
+                    m2 = m2[..., ::-1]
+                    # here we add dummy levels, might want to reconsider converting
+                    # "trimData" to "reOrderData" and use actual levels?
+                    m3 = numpy.concatenate(
+                        (m2, numpy.zeros(
+                            (m2.shape[0], 1))), axis=1)
+                except Exception:
+                    print("SETTING M3 to False")
+                    m3 = False  # no mesh but gengrid
+                    cellData = False
+
         # Could still be meshfill with mesh data
         # Ok probably should do a test for hgrid before sending data2
         if isinstance(gm, meshfill.Gfm) and data2 is not None:
@@ -382,15 +388,25 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False,
     if m3 is not None:
         # Create unstructured grid points
         vg = vtk.vtkUnstructuredGrid()
-        for i in range(numberOfCells):
-            pt_ids = []
-            for j in range(nVertices):
-                indx = i * nVertices + j
-                if not numpy.isnan(m3[indx][0]):  # missing value means skip vertex
-                    pt_ids.append(indx)
-            vg.InsertNextCell(vtk.VTK_POLYGON,
-                              len(pt_ids),
-                              pt_ids)
+        print("WELL M3 is:",m3)
+        if m3 is not False:
+            for i in range(numberOfCells):
+                pt_ids = []
+                for j in range(nVertices):
+                    indx = i * nVertices + j
+                    if not numpy.isnan(m3[indx][0]):  # missing value means skip vertex
+                        pt_ids.append(indx)
+                vg.InsertNextCell(vtk.VTK_POLYGON,
+                                  len(pt_ids),
+                                  pt_ids)
+        else:  # gengrid with no bounds
+            lat = g.getLatitude().asma()
+            lon = g.getLongitude().asma()
+            numberOfPoints = len(lat)
+            pts = vtk.vtkPoints()
+            for i in range(numberOfPoints):
+                pts.InsertNextPoint((lon[i],lat[i],0.))
+            vg.SetPoints(pts)
     else:
         # Ok a simple structured grid is enough
         if grid is None:
@@ -482,12 +498,14 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False,
 
     # attribute data
     gridForAttribute = grid if grid else vg
+    print("GFA AND VG:",gridForAttribute, vg)
     if genVectors:
         attribute = generateVectorArray(data1, data2, gridForAttribute)
     else:
         attribute = numpy_to_vtk_wrapper(data1.filled(0.).flat,
                                          deep=False)
         attribute.SetName("scalar")
+    print("CELLDATA:????",cellData)
     if cellData:
         attributes = gridForAttribute.GetCellData()
     else:
@@ -499,33 +517,38 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False,
 
     if grid is None:
         # First create the points/vertices (in vcs terms)
-        pts = vtk.vtkPoints()
-        # Convert nupmy array to vtk ones
-        ppV = numpy_to_vtk_wrapper(m3, deep=deep)
-        pts.SetData(ppV)
-        ptsBounds = pts.GetBounds()
-        xRange = ptsBounds[1] - ptsBounds[0]
-        xm, xM, ym, yM, tmp, tmp2 = pts.GetBounds()
-
-        # 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))
+        if m3 is not False:
+            pts = vtk.vtkPoints()
+            # Convert nupmy array to vtk ones
+            ppV = numpy_to_vtk_wrapper(m3, deep=deep)
+            pts.SetData(ppV)
+            ptsBounds = pts.GetBounds()
+            xRange = ptsBounds[1] - ptsBounds[0]
+            xm, xM, ym, yM, tmp, tmp2 = pts.GetBounds()
+
+            # 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))
 
-        vg.SetPoints(pts)
+            vg.SetPoints(pts)
+        else:
+            wc = [lon.min(),lon.max(),lat.min(),lat.max()]
+            print("WC:",wc)
         # index into the scalar array. Used for upgrading
         # the scalar after wrapping. Note this will work
         # correctly only for cell data. For point data
         # the indexes for points on the border will be incorrect after
         # wrapping
         pedigreeId = vtk.vtkIntArray()
         pedigreeId.SetName("PedigreeIds")
+        print("TUPLES:",attribute.GetNumberOfTuples())
         pedigreeId.SetNumberOfTuples(attribute.GetNumberOfTuples())
         for i in range(0, attribute.GetNumberOfTuples()):
             pedigreeId.SetValue(i, i)
@@ -539,11 +562,19 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False,
             vg = wrapDataSetX(vg)
             pts = vg.GetPoints()
             xm, xM, ym, yM, tmp, tmp2 = vg.GetPoints().GetBounds()
-        vg = doWrapData(vg, wc)
-        pts = vg.GetPoints()
         xm, xM, ym, yM, tmp, tmp2 = vg.GetPoints().GetBounds()
+        print("XS:",xm, xM, ym, yM, tmp, tmp2)
+        pts = vg.GetPoints()
+        print("HERE WE HAVE:",pts.GetNumberOfPoints())
+        if m3 is not False:
+            vg = doWrapData(vg, wc)
+            pts = vg.GetPoints()
+            xm, xM, ym, yM, tmp, tmp2 = vg.GetPoints().GetBounds()
         projection = vcs.elements["projection"][gm.projection]
-        vg.SetPoints(pts)
+        #vg.SetPoints(pts)
+        print("XS:",xm, xM, ym, yM, tmp, tmp2)
+        import pdb
+        pdb.set_trace()
         geo, geopts = project(pts, projection, getWrappedBounds(
             wc, [xm, xM, ym, yM], wrap))
         # proj4 returns inf for points that are not visible. Set those to a valid point