diff --git a/CMake/cmake_modules/FindMatplotlib.cmake b/CMake/cmake_modules/FindMatplotlib.cmake index 5cae4f77e0..fed7694e5d 100644 --- a/CMake/cmake_modules/FindMatplotlib.cmake +++ b/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 ) diff --git a/Packages/vcs/Lib/VTKPlots.py b/Packages/vcs/Lib/VTKPlots.py index b0afb987d9..e9150a2002 100644 --- a/Packages/vcs/Lib/VTKPlots.py +++ b/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]) diff --git a/Packages/vcs/Lib/vcs2vtk.py b/Packages/vcs/Lib/vcs2vtk.py index f237678061..be64b284d3 100644 --- a/Packages/vcs/Lib/vcs2vtk.py +++ b/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,6 +213,7 @@ 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: @@ -223,7 +221,7 @@ def genGridOnPoints(data1, gm, deep=True, grid=None, geo=None): 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 diff --git a/Packages/vcs/Lib/vcsvtk/vectorpipeline.py b/Packages/vcs/Lib/vcsvtk/vectorpipeline.py index 46f710fc59..f308522051 100644 --- a/Packages/vcs/Lib/vcsvtk/vectorpipeline.py +++ b/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) diff --git a/testing/vcs/CMakeLists.txt b/testing/vcs/CMakeLists.txt index 72bddad390..d4eb3cdb72 100644 --- a/testing/vcs/CMakeLists.txt +++ b/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() diff --git a/testing/vcs/test_vcs_vectors_robinson.py b/testing/vcs/test_vcs_vectors_robinson.py new file mode 100644 index 0000000000..1614fc03df --- /dev/null +++ b/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)