diff --git a/Packages/vcs/vcs/VTKPlots.py b/Packages/vcs/vcs/VTKPlots.py index 57e3f1057a..5a22482a37 100644 --- a/Packages/vcs/vcs/VTKPlots.py +++ b/Packages/vcs/vcs/VTKPlots.py @@ -755,11 +755,11 @@ def plotContinents(self, wc, projection, wrap, vp, priority, **kargs): if continents_path is None: return (None, 1, 1) contData = vcs2vtk.prepContinents(continents_path) + contData = vcs2vtk.doWrapData(contData, wc, fastClip=False) contMapper = vtk.vtkPolyDataMapper() contMapper.SetInputData(contData) contActor = vtk.vtkActor() contActor.SetMapper(contMapper) - contActor = vcs2vtk.doWrap(contActor, wc, fastClip=False) if projection.type != "linear": contData = contActor.GetMapper().GetInput() @@ -1471,24 +1471,19 @@ def update_input(self, vtkobjects, array1, array2=None, update=True): ports = vtkobjects["vtk_backend_geofilters"] else: # Vector plot + # TODO: this does not work with wrapping ports = vtkobjects["vtk_backend_glyphfilters"] w = vcs2vtk.generateVectorArray(array1, array2, vg) vg.GetPointData().AddArray(w) ports[0].SetInputData(vg) - projection = None - if "vtk_backend_geo" in vtkobjects: - projection = vtkobjects["vtk_backend_geo"] - if "vtk_backend_actors" in vtkobjects: i = 0 for a in vtkobjects["vtk_backend_actors"]: act = a[0] - wrp = a[1] if a[1] is missingMapper: i -= 1 mapper = missingMapper2 - wrp = a[2] else: # Labeled contours are a different kind if "vtk_backend_luts" in vtkobjects: @@ -1526,9 +1521,6 @@ def update_input(self, vtkobjects, array1, array2=None, update=True): mapper.SetScalarModeToUseCellData() mapper.SetScalarRange(rg[0], rg[1]) act.SetMapper(mapper) - if not projection: - # Wrap only if the data is not projected - act = vcs2vtk.doWrap(a[0], wrp) i += 1 taxis = array1.getTime() diff --git a/Packages/vcs/vcs/utils.py b/Packages/vcs/vcs/utils.py index 987db74754..c025b93ce2 100644 --- a/Packages/vcs/vcs/utils.py +++ b/Packages/vcs/vcs/utils.py @@ -1757,10 +1757,7 @@ def getworldcoordinates(gm, X, Y): return wc if (((not isinstance(X, cdms2.axis.TransientAxis) and isinstance(Y, cdms2.axis.TransientAxis)) or - not vcs.utils.monotonic(X[:])) and - numpy.allclose([datawc[0], datawc[1]], 1.e20))\ - or (hasattr(gm, "projection") and - vcs.elements["projection"][gm.projection].type != "linear"): + not vcs.utils.monotonic(X[:])) and numpy.allclose([datawc[0], datawc[1]], 1.e20)): wc[0] = X[:].min() wc[1] = X[:].max() if numpy.isclose(datawc[2], 1.e20): diff --git a/Packages/vcs/vcs/vcs2vtk.py b/Packages/vcs/vcs/vcs2vtk.py index 721c91f00e..e1ffbe070b 100644 --- a/Packages/vcs/vcs/vcs2vtk.py +++ b/Packages/vcs/vcs/vcs2vtk.py @@ -99,12 +99,11 @@ def putMaskOnVTKGrid(data, grid, actorColor=None, cellData=True, deep=True): if msk is not numpy.ma.nomask and not numpy.allclose(msk, False): if actorColor is not None: flatIMask = msk.astype(numpy.double).flat + grid2 = grid.NewInstance() if grid.IsA("vtkStructuredGrid"): - grid2 = vtk.vtkStructuredGrid() vtkmask = numpy_to_vtk_wrapper(flatIMask, deep=deep, array_type=vtk.VTK_DOUBLE) attributes2 = grid2.GetCellData() if cellData else grid2.GetPointData() else: - grid2 = vtk.vtkUnstructuredGrid() if (cellData): attributes2 = grid2.GetCellData() attributes = grid.GetCellData() @@ -134,7 +133,7 @@ def putMaskOnVTKGrid(data, grid, actorColor=None, cellData=True, deep=True): pointToCell.SetInputConnection(geoFilter.GetOutputPort()) geoFilter = pointToCell lut.SetNumberOfTableValues(256) - lut.SetTableValue(0, 1., 1., 1., 1.) + lut.SetTableValue(0, 1., 1., 1., 0.) for i in range(1, 256): lut.SetTableValue(i, r / 100., g / 100., b / 100., a / 100.) else: @@ -167,6 +166,12 @@ def putMaskOnVTKGrid(data, grid, actorColor=None, cellData=True, deep=True): attributes.AddArray(vtkghost) setArray(grid, ghost, vtk.vtkDataSetAttributes.GhostArrayName(), cellData, isScalars=False) + if (grid.GetExtentType() == vtk.VTK_PIECES_EXTENT): + if (cellData): + pass + else: + removeHiddenPoints(grid) + return mapper @@ -213,7 +218,7 @@ def getBoundsList(axis, hasCellData, dualGrid): def setInfToValid(geoPoints, ghost): ''' Set infinity points to a point that already exists in the list. - If a ghost array is passed, we also hide infinity points. + We also hide infinity points in the ghost array. We return true if any points are infinity ''' anyInfinity = False @@ -238,6 +243,32 @@ def setInfToValid(geoPoints, ghost): return anyInfinity +def removeHiddenPoints(grid): + ghost = grid.GetPointGhostArray() + if (not ghost): + return + pts = grid.GetPoints() + for i in range(pts.GetNumberOfPoints()): + if (ghost.GetValue(i) & vtk.vtkDataSetAttributes.HIDDENPOINT): + cells = vtk.vtkIdList() + # point hidden, remove all cells used by this point + grid.GetPointCells(i, cells) + for j in range(cells.GetNumberOfIds()): + grid.DeleteCell(cells.GetId(j)) + grid.RemoveDeletedCells() + + +def removeHiddenCells(grid): + ghost = grid.GetCellGhostArray() + if (not ghost): + return + for i in range(grid.GetNumberOfCells()): + if (ghost.GetValue(i) & vtk.vtkDataSetAttributes.HIDDENCELL): + # cell hidden, remove it + grid.DeleteCell(i) + grid.RemoveDeletedCells() + + def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, dualGrid=False): continents = False @@ -421,31 +452,7 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, ptsBounds = pts.GetBounds() xRange = ptsBounds[1] - ptsBounds[0] xm, xM, ym, yM, tmp, tmp2 = pts.GetBounds() - if (isinstance(g, cdms2.hgrid.TransientCurveGrid) and - xRange > 360 and not numpy.isclose(xRange, 360)): - 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 - # the indexes for points on the border will be incorrect after - # wrapping - pedigreeId = vtk.vtkIntArray() - pedigreeId.SetName("PedigreeIds") - pedigreeId.SetNumberOfTuples(attribute.GetNumberOfTuples()) - for i in range(0, attribute.GetNumberOfTuples()): - pedigreeId.SetValue(i, i) - if cellData: - vg.GetCellData().SetPedigreeIds(pedigreeId) - else: - vg.GetPointData().SetPedigreeIds(pedigreeId) - vg = wrapDataSetX(vg) - pts = vg.GetPoints() - xm, xM, ym, yM, tmp, tmp2 = vg.GetPoints().GetBounds() - else: - xm, xM, ym, yM, tmp, tmp2 = grid.GetPoints().GetBounds() - projection = vcs.elements["projection"][gm.projection] - if grid is None: - vg.SetPoints(pts) + # 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 @@ -456,6 +463,33 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, 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 + # the indexes for points on the border will be incorrect after + # wrapping + pedigreeId = vtk.vtkIntArray() + pedigreeId.SetName("PedigreeIds") + pedigreeId.SetNumberOfTuples(attribute.GetNumberOfTuples()) + for i in range(0, attribute.GetNumberOfTuples()): + pedigreeId.SetValue(i, i) + if cellData: + vg.GetCellData().SetPedigreeIds(pedigreeId) + else: + vg.GetPointData().SetPedigreeIds(pedigreeId) + + if (isinstance(g, cdms2.hgrid.TransientCurveGrid) and + xRange > 360 and not numpy.isclose(xRange, 360)): + 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() + projection = vcs.elements["projection"][gm.projection] + 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 @@ -477,9 +511,15 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, ym = p[1] if (p[1] > yM): yM = p[1] + # hidden point don't work for polys or unstructured grids. + # We remove the cells in this case. + if (vg.GetExtentType() == vtk.VTK_PIECES_EXTENT): + removeHiddenPoints(vg) + # Sets the vertics into the grid vg.SetPoints(geopts) else: + xm, xM, ym, yM, tmp, tmp2 = grid.GetPoints().GetBounds() vg = grid # Add a GlobalIds array to keep track of cell ids throughout the pipeline globalIds = numpy_to_vtk_wrapper(numpy.arange(0, vg.GetNumberOfCells()), deep=True) @@ -855,32 +895,46 @@ def dump2VTK(obj, fnm=None): dsw.Write() -# Wrapping around -def doWrap(Act, wc, wrap=[0., 360], fastClip=True): +def doWrapData(data, wc, wrap=[0., 360], fastClip=True): + ''' + Wrapping around and 'wrap' modulo' and clipping. + wrap contains YWrap, XWrap modulo in degrees, 0 means no wrap + ''' if wrap is None: - return Act - Mapper = Act.GetMapper() - Mapper.Update() - data = Mapper.GetInput() + return data + + # convert to poly data + surface = vtk.vtkDataSetSurfaceFilter() + surface.SetInputData(data) + surface.Update() + data = surface.GetOutput() + + bounds = data.GetBounds() # insure that GLOBALIDS are not removed by the append filter attributes = data.GetCellData() attributes.SetActiveAttribute(-1, attributes.GLOBALIDS) xmn = min(wc[0], wc[1]) xmx = max(wc[0], wc[1]) - if numpy.allclose(xmn, 1.e20) or numpy.allclose(xmx, 1.e20): - xmx = abs(wrap[1]) - xmn = -wrap[1] + if (numpy.allclose(xmn, 1.e20) or numpy.allclose(xmx, 1.e20)): + xmn = bounds[0] + if (wrap[1] > 0. and (bounds[1] - bounds[0] > wrap[1])): + xmx = bounds[0] + abs(wrap[1]) + else: + xmx = bounds[1] ymn = min(wc[2], wc[3]) ymx = max(wc[2], wc[3]) if numpy.allclose(ymn, 1.e20) or numpy.allclose(ymx, 1.e20): - ymx = abs(wrap[0]) - ymn = -wrap[0] + ymn = bounds[2] + if (wrap[0] > 0. and (bounds[3] - bounds[2] > wrap[0])): + ymx = bounds[2] + abs(wrap[0]) + else: + ymx = bounds[3] appendFilter = vtk.vtkAppendPolyData() appendFilter.AddInputData(data) appendFilter.Update() # X axis wrappping - Amn, Amx = Act.GetXRange() + Amn, Amx = bounds[0], bounds[1] if wrap[1] != 0.: i = 0 while Amn > xmn: @@ -908,7 +962,7 @@ def doWrap(Act, wc, wrap=[0., 360], fastClip=True): appendFilter.Update() # Y axis wrapping - Amn, Amx = Act.GetYRange() + Amn, Amx = bounds[2], bounds[3] if wrap[0] != 0.: i = 0 while Amn > ymn: @@ -957,8 +1011,7 @@ def doWrap(Act, wc, wrap=[0., 360], fastClip=True): attributes.GetArray("GlobalIds", globalIdsIndex) attributes.SetActiveAttribute(globalIdsIndex, attributes.GLOBALIDS) - Mapper.SetInputData(clipper.GetOutput()) - return Act + return clipper.GetOutput() # Wrap grid in interval minX, minX + 360 diff --git a/Packages/vcs/vcs/vcsvtk/boxfillpipeline.py b/Packages/vcs/vcs/vcsvtk/boxfillpipeline.py index 005241b4a0..0847b318b4 100644 --- a/Packages/vcs/vcs/vcsvtk/boxfillpipeline.py +++ b/Packages/vcs/vcs/vcsvtk/boxfillpipeline.py @@ -1,5 +1,4 @@ from .pipeline2d import Pipeline2D -from .. import vcs2vtk import fillareautils import numpy @@ -84,12 +83,6 @@ def _plotInternal(self): act = vtk.vtkActor() act.SetMapper(mapper) - if self._vtkGeoTransform is None: - # If using geofilter on wireframed does not get wrppaed not - # sure why so sticking to many mappers - act = vcs2vtk.doWrap(act, plotting_dataset_bounds, - self._dataWrapModulo) - # TODO We shouldn't need this conditional branch, the 'else' body # should be used and GetMapper called to get the mapper as needed. # If this is needed for other reasons, we need a comment explaining diff --git a/Packages/vcs/vcs/vcsvtk/isofillpipeline.py b/Packages/vcs/vcs/vcsvtk/isofillpipeline.py index 273376c090..c8c42f192b 100644 --- a/Packages/vcs/vcs/vcsvtk/isofillpipeline.py +++ b/Packages/vcs/vcs/vcsvtk/isofillpipeline.py @@ -1,5 +1,4 @@ from .pipeline2d import Pipeline2D -from .. import vcs2vtk import fillareautils import numpy @@ -124,12 +123,6 @@ def _plotInternal(self): act = vtk.vtkActor() act.SetMapper(mapper) - if self._vtkGeoTransform is None: - # If using geofilter on wireframed does not get wrppaed not - # sure why so sticking to many mappers - act = vcs2vtk.doWrap(act, plotting_dataset_bounds, - self._dataWrapModulo) - patact = None # TODO see comment in boxfill. if mapper is self._maskedDataMapper: @@ -163,7 +156,7 @@ def _plotInternal(self): wc=plotting_dataset_bounds, geoBounds=self._vtkDataSet.GetBounds(), geo=self._vtkGeoTransform, priority=self._template.data.priority, - create_renderer=(dataset_renderer is None)) + create_renderer=(mapper is self._maskedDataMapper or dataset_renderer is None)) for act in patternActors: self._context().fitToViewport( act, vp, diff --git a/Packages/vcs/vcs/vcsvtk/isolinepipeline.py b/Packages/vcs/vcs/vcsvtk/isolinepipeline.py index 02a55c4cbc..bb4f7e3c58 100644 --- a/Packages/vcs/vcs/vcsvtk/isolinepipeline.py +++ b/Packages/vcs/vcs/vcsvtk/isolinepipeline.py @@ -251,11 +251,6 @@ def _plotInternal(self): p.SetLineWidth(tmpLineWidths[i]) vcs2vtk.stippleLine(p, tmpLineTypes[i]) - if self._vtkGeoTransform is None: - # If using geofilter on wireframed does not get wrppaed not - # sure why so sticking to many mappers - act = vcs2vtk.doWrap(act, plotting_dataset_bounds, - self._dataWrapModulo) actors.append([act, plotting_dataset_bounds]) # create a new renderer for this mapper @@ -283,11 +278,6 @@ def _plotInternal(self): mappers.insert(0, self._maskedDataMapper) act = vtk.vtkActor() act.SetMapper(self._maskedDataMapper) - if self._vtkGeoTransform is None: - # If using geofilter on wireframed does not get wrppaed not - # sure why so sticking to many mappers - act = vcs2vtk.doWrap(act, plotting_dataset_bounds, - self._dataWrapModulo) actors.append([act, self._maskedDataMapper, plotting_dataset_bounds]) # create a new renderer for this mapper # (we need one for each mapper because of cmaera flips) diff --git a/Packages/vcs/vcs/vcsvtk/meshfillpipeline.py b/Packages/vcs/vcs/vcsvtk/meshfillpipeline.py index 3bf2bb72d1..4f28105f5f 100644 --- a/Packages/vcs/vcs/vcsvtk/meshfillpipeline.py +++ b/Packages/vcs/vcs/vcsvtk/meshfillpipeline.py @@ -175,12 +175,6 @@ def _plotInternal(self): prop = act.GetProperty() prop.SetRepresentationToWireframe() - if self._vtkGeoTransform is None: - # If using geofilter on wireframed does not get wrppaed not - # sure why so sticking to many mappers - act = vcs2vtk.doWrap(act, plotting_dataset_bounds, - self._dataWrapModulo) - # TODO See comment in boxfill. if mapper is self._maskedDataMapper: actors.append([act, self._maskedDataMapper, plotting_dataset_bounds]) diff --git a/testing/vcs/CMakeLists.txt b/testing/vcs/CMakeLists.txt index fe9cfa412f..05b2f1ebf6 100644 --- a/testing/vcs/CMakeLists.txt +++ b/testing/vcs/CMakeLists.txt @@ -14,12 +14,25 @@ cdat_add_test(test_vcs_bad_png_path ) foreach(projection polar mollweide lambert orthographic mercator polyconic robinson) - cdat_add_test(test_vcs_boxfill_${projection} - "${PYTHON_EXECUTABLE}" - ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_projection.py - "${BASELINE_DIR}/test_vcs_boxfill_${projection}.png" - ${projection} - ) + foreach(zoom none datawc subset) + if (zoom STREQUAL "none") + cdat_add_test(test_vcs_boxfill_${projection} + "${PYTHON_EXECUTABLE}" + ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_projection.py + "${BASELINE_DIR}/test_vcs_boxfill_${projection}.png" + ${projection} + ${zoom} + ) + else() + cdat_add_test(test_vcs_boxfill_${projection}_${zoom} + "${PYTHON_EXECUTABLE}" + ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_projection.py + "${BASELINE_DIR}/test_vcs_boxfill_${projection}_${zoom}.png" + ${projection} + ${zoom} + ) + endif() + endforeach() endforeach() foreach(lat_0 45 90) @@ -783,11 +796,15 @@ cdat_add_test(test_vcs_settings_color_name_rgba # ${cdat_SOURCE_DIR}/testing/vcs/test_aspect_ratio.py # ${cdat_SOURCE_DIR}/testing/vcs/test_aspect_ratio.py # ) - cdat_add_test(test_vcs_polar_set_opt_param_polar - "${PYTHON_EXECUTABLE}" - ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_polar_set_opt_param_polar.py - "${BASELINE_DIR}/test_vcs_polar_set_opt_param_polar.png" - ) + + foreach(z none subset datawc0 datawc1 datawc2) + cdat_add_test(test_vcs_polar_zoom_${z} + "${PYTHON_EXECUTABLE}" + ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_polar_zoom.py + "${BASELINE_DIR}/test_vcs_polar_zoom_${z}.png" + ${z}) + endforeach() + cdat_add_test(test_vcs_boxfill_lev1_lev2 "${PYTHON_EXECUTABLE}" ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_lev1_lev2.py diff --git a/testing/vcs/test_vcs_boxfill_projection.py b/testing/vcs/test_vcs_boxfill_projection.py index 6f319efd4b..1f262013bf 100644 --- a/testing/vcs/test_vcs_boxfill_projection.py +++ b/testing/vcs/test_vcs_boxfill_projection.py @@ -2,7 +2,7 @@ baselineName = sys.argv[1] projection = sys.argv[2] - +zoom = sys.argv[3] f = cdms2.open(vcs.sample_data + "/clt.nc") a = f("clt") @@ -11,7 +11,16 @@ p = x.getprojection(projection) b = x.createboxfill() b.projection = p -x.plot(a(latitude=(90,-90)), b, bg=1) +if (zoom == 'none'): + x.plot(a(latitude=(90,-90)), b, bg=1) +elif (zoom == 'subset'): + x.plot(a(latitude=(-50,90), longitude=(30, -30)), b, bg=1) +else: + b.datawc_x1 = 30 + b.datawc_x2 = -30 + b.datawc_y1 = -50 + b.datawc_y2 = 90 + x.plot(a, b, bg=1) fileName = os.path.basename(baselineName) fileName = os.path.splitext(fileName)[0] diff --git a/testing/vcs/test_vcs_meshfill_draw_mesh.py b/testing/vcs/test_vcs_meshfill_draw_mesh.py index 08801d7a6d..a048811a96 100644 --- a/testing/vcs/test_vcs_meshfill_draw_mesh.py +++ b/testing/vcs/test_vcs_meshfill_draw_mesh.py @@ -8,4 +8,4 @@ m.mesh = True x.plot(s,m,bg=1) -regression.run(x, "test_meshfill_draw_mesh.png") \ No newline at end of file +regression.run(x, "test_vcs_meshfill_draw_mesh.png") diff --git a/testing/vcs/test_vcs_polar_set_opt_param_polar.py b/testing/vcs/test_vcs_polar_set_opt_param_polar.py deleted file mode 100644 index 4e777fb2b3..0000000000 --- a/testing/vcs/test_vcs_polar_set_opt_param_polar.py +++ /dev/null @@ -1,13 +0,0 @@ -import vcs, cdms2, sys, os, testing.regression as regression - -f = cdms2.open(os.path.join(vcs.sample_data,'clt.nc')) -s = f("clt",slice(0,1),squeeze=1) -x = regression.init() -x.setantialiasing(0) -x.drawlogooff() -x.setbgoutputdimensions(1200,1091,units="pixels") -i = x.createisofill() -p = x.getprojection("polar") -i.projection=p -x.plot(s,i,bg=1) -regression.run(x, "test_polar_set_opt_param_polar.png") \ No newline at end of file diff --git a/testing/vcs/test_vcs_polar_zoom.py b/testing/vcs/test_vcs_polar_zoom.py new file mode 100644 index 0000000000..5e0326c423 --- /dev/null +++ b/testing/vcs/test_vcs_polar_zoom.py @@ -0,0 +1,27 @@ +import vcs, cdms2, sys, os, testing.regression as regression + +zoom = sys.argv[2] + +f = cdms2.open(os.path.join(vcs.sample_data,'clt.nc')) +s = f("clt",slice(0,1),squeeze=1) +x = regression.init() +i = x.createisofill() +p = x.getprojection("polar") +i.projection=p +if (zoom == 'none'): + x.plot(s,i,bg=1) +elif (zoom == 'subset'): + x.plot(s(latitude=(-50,90), longitude=(30, -30)), i, bg=1) +else: + i.datawc_x1 = 30 + i.datawc_x2 = -30 + i.datawc_y1 = -50 + i.datawc_y2 = 90 + if (zoom == 'datawc1'): + i.datawc_x1, i.datawc_x2 = i.datawc_x2, i.datawc_x1 + if (zoom == 'datawc2'): + i.datawc_y1, i.datawc_y2 = i.datawc_y2, i.datawc_y1 + x.plot(s, i, bg=1) + +file = "test_vcs_polar_zoom_" + zoom + ".png" +regression.run(x, file)