From 2cfa99adc592243871514ed2e704c01c14b9f507 Mon Sep 17 00:00:00 2001 From: Charles Doutriaux Date: Wed, 28 Feb 2018 16:05:15 -0800 Subject: [PATCH 1/2] painfully trying to plot isofill with gengrid --- vcs/Canvas.py | 4 +- vcs/vcs2vtk.py | 123 +++++++++++++++++++++++++++++++------------------ 2 files changed, 79 insertions(+), 48 deletions(-) diff --git a/vcs/Canvas.py b/vcs/Canvas.py index 5d451c351..81be57b12 100644 --- a/vcs/Canvas.py +++ b/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' % diff --git a/vcs/vcs2vtk.py b/vcs/vcs2vtk.py index 9fc53e064..7bf0a77a7 100644 --- a/vcs/vcs2vtk.py +++ b/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,26 +517,30 @@ 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 @@ -526,6 +548,7 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, # 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 From bd6349d95f64eec22e98aeba8cd630c7e36b8ac1 Mon Sep 17 00:00:00 2001 From: Charles Doutriaux Date: Mon, 9 Apr 2018 10:54:59 -0700 Subject: [PATCH 2/2] blah --- vcs/VTKPlots.py | 9 ++- vcs/vcs2vtk.py | 116 +++++++++++++++++++-------------------- vcs/vcsvtk/pipeline2d.py | 15 ++++- 3 files changed, 75 insertions(+), 65 deletions(-) diff --git a/vcs/VTKPlots.py b/vcs/VTKPlots.py index fa5c55d30..00711bb99 100644 --- a/vcs/VTKPlots.py +++ b/vcs/VTKPlots.py @@ -1027,10 +1027,13 @@ def trimData2D(self, data): g = data.getGrid() gaxes = list(g.getAxisList()) daxes = list(data.getAxisList()) - if daxes[len(daxes) - len(gaxes):] == gaxes: + if numpy.allclose(daxes[len(daxes) - len(gaxes):], gaxes): # Ok it is gridded and the grid axes are last - return self.cleanupData( - data(*(slice(0, 1),) * (len(daxes) - len(gaxes)), squeeze=1)) + if len(daxes) != len(gaxes): + return self.cleanupData( + data(*(slice(0, 1),) * (len(daxes) - len(gaxes)), squeeze=1)) + else: + return data else: # Ok just return the last two dims return self.cleanupData( diff --git a/vcs/vcs2vtk.py b/vcs/vcs2vtk.py index dae19fc3e..4a672f271 100644 --- a/vcs/vcs2vtk.py +++ b/vcs/vcs2vtk.py @@ -356,10 +356,33 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, 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 + except Exception: # ok no mesh, so unstructured grid point data + vg = vtk.vtkUnstructuredGrid() + lat = g.getLatitude().asma() + lon = g.getLongitude().asma() + pts = vtk.vtkPoints() + for i in range(len(lat)): + pts.InsertNextPoint(float(lon[i]),float(lat[i]),0.) + vg.SetPoints(pts) + + + + xm, xM, ym, yM = lon.min(),lon.max(),lat.min(),lat.max() + wc = [gm.datawc_x1, gm.datawc_x2, gm.datawc_y1, gm.datawc_y2] + geo, geopts = project(pts, projection, getWrappedBounds(wc, [xm, xM, ym, yM], wrap)) + out = {"vtk_backend_grid": vg, + "xm": xm, + "xM": xM, + "ym": ym, + "yM": yM, + "continents": continents, + "wrap": wrap, + "geo": geo, + "cellData": False, + "data": data1, + "data2": data2 + } + return out # Could still be meshfill with mesh data # Ok probably should do a test for hgrid before sending data2 @@ -388,25 +411,15 @@ 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() - 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) + 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: # Ok a simple structured grid is enough if grid is None: @@ -498,14 +511,12 @@ 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: @@ -517,30 +528,26 @@ 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) - 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) + 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 = [lon.min(),lon.max(),lat.min(),lat.max()] - print("WC:",wc) + wc = vcs.utils.getworldcoordinates(gm, + data1.getAxis(-1), + data1.getAxis(-2)) + + vg.SetPoints(pts) # index into the scalar array. Used for upgrading # the scalar after wrapping. Note this will work # correctly only for cell data. For point data @@ -548,7 +555,6 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, # 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) @@ -565,18 +571,8 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, vg = doWrapData(vg, wc, wrap) 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) - print("XS:",xm, xM, ym, yM, tmp, tmp2) - import pdb - pdb.set_trace() + vg.SetPoints(pts) 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 diff --git a/vcs/vcsvtk/pipeline2d.py b/vcs/vcsvtk/pipeline2d.py index 54f1a9ed1..37d06dc74 100644 --- a/vcs/vcsvtk/pipeline2d.py +++ b/vcs/vcsvtk/pipeline2d.py @@ -289,7 +289,7 @@ def plot(self, data1, data2, tmpl, grid, transform, **kargs): self._updateScalarData() self._min = self._data1.min() self._max = self._data1.max() - self._scalarRange = vcs.minmax(self._data1) + self._scalarRange = (self._min, self._max) # Create/update the VTK dataset. plotBasedDualGrid = kargs.get('plot_based_dual_grid', True) @@ -357,7 +357,18 @@ def _updateVTKDataSet(self, plotBasedDualGrid): def _createPolyDataFilter(self): """This is only used when we use the grid stored in the file for all plots.""" self._vtkPolyDataFilter = vtk.vtkDataSetSurfaceFilter() - if self._hasCellData == self._needsCellData: + import ipdb + ipdb.set_trace() + if self._vtkDataSet.IsA("vtkUnstructuredGrid"): + print("UNSADFRRTVXDFGFDGDHRFDCGDRFDSGDSFGREDSCFDGRDFSG") + delny = vtk.vtkDelaunay2D() + profile = vtk.vtkPolyData() + profile.SetPoints(self._vtkDataSet.GetPoints()) + delny.SetInputData(profile) + delny.SetTolerance(0.001) + self._vtkPolyDataFilter.SetInputConnection(delny.GetOutputPort()) + + elif self._hasCellData == self._needsCellData: self._vtkPolyDataFilter.SetInputData(self._vtkDataSet) elif self._hasCellData: # use cells but needs points