diff --git a/Packages/vcs/vcs/Canvas.py b/Packages/vcs/vcs/Canvas.py index 49361248ce..e95609e3ea 100644 --- a/Packages/vcs/vcs/Canvas.py +++ b/Packages/vcs/vcs/Canvas.py @@ -71,7 +71,7 @@ import vcsaddons # noqa import vcs.manageElements # noqa import configurator # noqa -from projection import round_projections # noqa +from projection import no_deformation_projections # noqa # Python < 3 DeprecationWarning ignored by default warnings.simplefilter('default') @@ -3497,7 +3497,7 @@ def set_convert_labels(copy_mthd, test=0): if hasattr(gm, "priority") and gm.priority == 0: return p = self.getprojection(gm.projection) - if p.type in round_projections and ( + if p.type in no_deformation_projections and ( doratio == "0" or doratio[:4] == "auto"): doratio = "1t" for keyarg in keyargs.keys(): @@ -3554,7 +3554,7 @@ def set_convert_labels(copy_mthd, test=0): t.data.y2 = p.viewport[3] proj = self.getprojection(p.projection) - if proj.type in round_projections and ( + if proj.type in no_deformation_projections and ( doratio == "0" or doratio[:4] == "auto"): doratio = "1t" @@ -3610,7 +3610,7 @@ def set_convert_labels(copy_mthd, test=0): tp = "textcombined" gm = vcs.elements[tp][arglist[4]] p = self.getprojection(gm.projection) - if p.type in round_projections: + if p.type in no_deformation_projections: doratio = "1t" if p.type == 'linear': if gm.g_name == 'Gfm': diff --git a/Packages/vcs/vcs/VTKPlots.py b/Packages/vcs/vcs/VTKPlots.py index 814719536b..9d3d85c748 100644 --- a/Packages/vcs/vcs/VTKPlots.py +++ b/Packages/vcs/vcs/VTKPlots.py @@ -597,6 +597,7 @@ def plot(self, data1, data2, template, gtype, gname, bg, *args, **kargs): vtk_backend_grid = kargs.get("vtk_backend_grid", None) vtk_backend_geo = kargs.get("vtk_backend_geo", None) + bounds = vtk_backend_grid.GetBounds() if vtk_backend_grid else None pipeline = vcsvtk.createPipeline(gm, self) if pipeline is not None: @@ -626,7 +627,7 @@ def plot(self, data1, data2, template, gtype, gname, bg, *args, **kargs): ren, to=to, tt=tt, - cmap=self.canvas.colormap) + cmap=self.canvas.colormap, geoBounds=bounds, geo=vtk_backend_geo) self.setLayer(ren, tt.priority) self.text_renderers[tt_key] = ren elif gtype == "line": @@ -635,7 +636,6 @@ def plot(self, data1, data2, template, gtype, gname, bg, *args, **kargs): cmap=self.canvas.colormap) returned["vtk_backend_line_actors"] = actors create_renderer = True - bounds = vtk_backend_grid.GetBounds() if vtk_backend_grid else None for act, geo in actors: ren = self.fitToViewport( act, diff --git a/Packages/vcs/vcs/projection.py b/Packages/vcs/vcs/projection.py index a8476a9892..6b19c5e62f 100644 --- a/Packages/vcs/vcs/projection.py +++ b/Packages/vcs/vcs/projection.py @@ -16,13 +16,15 @@ import vcs import copy -# projection that seems to be doing a circle -# We will probably to add some more in it as we find more that fit this -round_projections = ['polar (non gctp)', 'stereographic', - 'orthographic', "ortho", ] +# used to decide if we show longitude labels for round projections or +# latitude labels for elliptical projections +round_projections = ['polar (non gctp)', 'stereographic'] +elliptical_projections = ["robinson", "mollweide", 'orthographic', "ortho"] +# projections in this list are not deformed based on the window size +no_deformation_projections = ['polar (non gctp)', 'stereographic', + 'orthographic', "ortho", ] no_over_proj4_parameter_projections = round_projections+["aeqd", "lambert conformal c"] -elliptical_projections = ["robinson", "mollweide"] def process_src(nm, code): diff --git a/Packages/vcs/vcs/vcs2vtk.py b/Packages/vcs/vcs/vcs2vtk.py index 81142492c5..686b1a6774 100644 --- a/Packages/vcs/vcs/vcs2vtk.py +++ b/Packages/vcs/vcs/vcs2vtk.py @@ -4,12 +4,14 @@ import numpy import json import os +import math import meshfill from vtk.util import numpy_support as VN import cdms2 import warnings from projection import round_projections, no_over_proj4_parameter_projections from vcsvtk import fillareautils +import sys import numbers f = open(os.path.join(vcs.prefix, "share", "vcs", "wmo_symbols.json")) @@ -163,23 +165,6 @@ def putMaskOnVTKGrid(data, grid, actorColor=None, cellData=True, deep=True): return mapper -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. - ptype = projDict.get(projection._type, projection.type) - if (ptype.lower() == "merc"): - lat = data.getLatitude() - if isinstance(lat, cdms2.axis.TransientAxis): - lat = lat[:] - # Reverse the latitudes incase the starting latitude is greater - # than the ending one - if lat[-1] < lat[0]: - lat = lat[::-1] - data = data(latitude=(max(-85, lat.min()), min(85, lat.max()))) - return data - - def getBoundsList(axis, hasCellData, dualGrid): ''' Returns the bounds list for 'axis'. If axis has n elements the @@ -220,6 +205,34 @@ def getBoundsList(axis, hasCellData, dualGrid): return None +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 return true if any points are infinity + ''' + anyInfinity = False + validPoint = [0, 0, 0] + for i in range(geoPoints.GetNumberOfPoints()): + point = geoPoints.GetPoint(i) + if (not math.isinf(point[0]) and not math.isinf(point[1])): + validPoint[0] = point[0] + validPoint[1] = point[1] + break + for i in range(geoPoints.GetNumberOfPoints()): + point = geoPoints.GetPoint(i) + if (math.isinf(point[0]) or math.isinf(point[1])): + anyInfinity = True + newPoint = list(point) + if (math.isinf(point[0])): + newPoint[0] = validPoint[0] + if (math.isinf(point[1])): + newPoint[1] = validPoint[1] + geoPoints.SetPoint(i, newPoint) + ghost.SetValue(i, vtk.vtkDataSetAttributes.HIDDENPOINT) + return anyInfinity + + def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, dualGrid=False): continents = False @@ -230,10 +243,6 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, xm, xM, ym, yM = None, None, None, None projection = vcs.elements["projection"][gm.projection] - data1 = handleProjectionEdgeCases(projection, data1) - if data2 is not None: - data2 = handleProjectionEdgeCases(projection, data2) - try: # First try to see if we can get a mesh out of this g = data1.getGrid() # Ok need unstructured grid @@ -444,6 +453,25 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False, data1.getAxis(-2)) 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 + # and hide them. + ghost = vg.AllocatePointGhostArray() + if (setInfToValid(geopts, ghost)): + # if there are hidden points, we recompute the bounds + xm = ym = sys.float_info.max + xM = yM = - sys.float_info.max + for i in range(pts.GetNumberOfPoints()): + if (ghost.GetValue(i) & vtk.vtkDataSetAttributes.HIDDENPOINT == 0): + # point not hidden + p = pts.GetPoint(i) + if (p[0] < xm): + xm = p[0] + if (p[0] > xM): + xM = p[0] + if (p[1] < ym): + ym = p[1] + if (p[1] > yM): + yM = p[1] # Sets the vertics into the grid vg.SetPoints(geopts) else: @@ -557,24 +585,42 @@ def apply_proj_parameters(pd, projection, x1, x2, y1, y2): else: pd.SetOptionalParameter("over", "false") setProjectionParameters(pd, projection) - if (hasattr(projection, 'centralmeridian') and - numpy.allclose(projection.centralmeridian, 1e+20)): - pd.SetCentralMeridian(float(x1 + x2) / 2.0) - if (hasattr(projection, 'centerlongitude') and - numpy.allclose(projection.centerlongitude, 1e+20)): - pd.SetOptionalParameter("lon_0", str(float(x1 + x2) / 2.0)) - if (hasattr(projection, 'originlatitude') and - numpy.allclose(projection.originlatitude, 1e+20)): - pd.SetOptionalParameter("lat_0", str(float(y1 + y2) / 2.0)) - if (hasattr(projection, 'centerlatitude') and - numpy.allclose(projection.centerlatitude, 1e+20)): - pd.SetOptionalParameter("lat_0", str(float(y1 + y2) / 2.0)) - if (hasattr(projection, 'standardparallel1') and - numpy.allclose(projection.standardparallel1, 1.e20)): - pd.SetOptionalParameter('lat_1', str(min(y1, y2))) - if (hasattr(projection, 'standardparallel2') and - numpy.allclose(projection.standardparallel2, 1.e20)): - pd.SetOptionalParameter('lat_2', str(max(y1, y2))) + if (hasattr(projection, 'centralmeridian')): + if (numpy.allclose(projection.centralmeridian, 1e+20)): + centralmeridian = float(x1 + x2) / 2.0 + else: + centralmeridian = projection.centralmeridian + pd.SetCentralMeridian(centralmeridian) + if (hasattr(projection, 'centerlongitude')): + if (numpy.allclose(projection.centerlongitude, 1e+20)): + centerlongitude = float(x1 + x2) / 2.0 + else: + centerlongitude = projection.centerlongitude + pd.SetOptionalParameter("lon_0", str(centerlongitude)) + if (hasattr(projection, 'originlatitude')): + if (numpy.allclose(projection.originlatitude, 1e+20)): + originlatitude = float(y1 + y2) / 2.0 + else: + originlatitude = projection.originlatitude + pd.SetOptionalParameter("lat_0", str(originlatitude)) + if (hasattr(projection, 'centerlatitude')): + if (numpy.allclose(projection.centerlatitude, 1e+20)): + centerlatitude = float(y1 + y2) / 2.0 + else: + centerlatitude = projection.centerlatitude + pd.SetOptionalParameter("lat_0", str(centerlatitude)) + if (hasattr(projection, 'standardparallel1')): + if (numpy.allclose(projection.standardparallel1, 1.e20)): + standardparallel1 = min(y1, y2) + else: + standardparallel1 = projection.standardparallel1 + pd.SetOptionalParameter('lat_1', str(standardparallel1)) + if (hasattr(projection, 'standardparallel2')): + if (numpy.allclose(projection.standardparallel2, 1.e20)): + standardparallel2 = max(y1, y2) + else: + standardparallel2 = projection.standardparallel2 + pd.SetOptionalParameter('lat_2', str(standardparallel2)) def projectArray(w, projection, wc, geo=None): @@ -1072,7 +1118,7 @@ def prepTextProperty(p, winSize, to="default", tt="default", cmap=None, def genTextActor(renderer, string=None, x=None, y=None, - to='default', tt='default', cmap=None): + to='default', tt='default', cmap=None, geoBounds=None, geo=None): if isinstance(to, str): to = vcs.elements["textorientation"][to] if isinstance(tt, str): @@ -1096,21 +1142,8 @@ def genTextActor(renderer, string=None, x=None, y=None, sz = renderer.GetRenderWindow().GetSize() actors = [] pts = vtk.vtkPoints() - geo = None if vcs.elements["projection"][tt.projection].type != "linear": - # Need to figure out new WC - Npts = 20 - for i in range(Npts + 1): - X = tt.worldcoordinate[ - 0] + float(i) / Npts * (tt.worldcoordinate[1] - - tt.worldcoordinate[0]) - for j in range(Npts + 1): - Y = tt.worldcoordinate[ - 2] + float(j) / Npts * (tt.worldcoordinate[3] - - tt.worldcoordinate[2]) - pts.InsertNextPoint(X, Y, 0.) - geo, pts = project(pts, tt.projection, tt.worldcoordinate, geo=None) - wc = pts.GetBounds()[:4] + wc = geoBounds[:4] # renderer.SetViewport(tt.viewport[0],tt.viewport[2],tt.viewport[1],tt.viewport[3]) renderer.SetWorldPoint(wc) @@ -1120,8 +1153,8 @@ def genTextActor(renderer, string=None, x=None, y=None, prepTextProperty(p, sz, to, tt, cmap) pts = vtk.vtkPoints() pts.InsertNextPoint(x[i], y[i], 0.) - if geo is not None: - geo, pts = project(pts, tt.projection, tt.worldcoordinate, geo=geo) + if vcs.elements["projection"][tt.projection].type != "linear": + _, pts = project(pts, tt.projection, tt.worldcoordinate, geo=geo) X, Y, tz = pts.GetPoint(0) X, Y = world2Renderer(renderer, X, Y, tt.viewport, wc) else: diff --git a/Packages/vcs/vcs/vcsvtk/boxfillpipeline.py b/Packages/vcs/vcs/vcsvtk/boxfillpipeline.py index f2a3ea6020..005241b4a0 100644 --- a/Packages/vcs/vcs/vcsvtk/boxfillpipeline.py +++ b/Packages/vcs/vcs/vcsvtk/boxfillpipeline.py @@ -152,7 +152,8 @@ def _plotInternal(self): z = None kwargs = {"vtk_backend_grid": self._vtkDataSet, "dataset_bounds": self._vtkDataSetBounds, - "plotting_dataset_bounds": plotting_dataset_bounds} + "plotting_dataset_bounds": plotting_dataset_bounds, + "vtk_backend_geo": self._vtkGeoTransform} if ("ratio_autot_viewport" in self._resultDict): kwargs["ratio_autot_viewport"] = vp self._resultDict.update(self._context().renderTemplate( diff --git a/Packages/vcs/vcs/vcsvtk/isofillpipeline.py b/Packages/vcs/vcs/vcsvtk/isofillpipeline.py index 55098f9e5c..887c6158c9 100644 --- a/Packages/vcs/vcs/vcsvtk/isofillpipeline.py +++ b/Packages/vcs/vcs/vcsvtk/isofillpipeline.py @@ -176,7 +176,8 @@ def _plotInternal(self): z = None kwargs = {"vtk_backend_grid": self._vtkDataSet, "dataset_bounds": self._vtkDataSetBounds, - "plotting_dataset_bounds": plotting_dataset_bounds} + "plotting_dataset_bounds": plotting_dataset_bounds, + "vtk_backend_geo": self._vtkGeoTransform} if ("ratio_autot_viewport" in self._resultDict): kwargs["ratio_autot_viewport"] = vp self._resultDict.update(self._context().renderTemplate( diff --git a/Packages/vcs/vcs/vcsvtk/isolinepipeline.py b/Packages/vcs/vcs/vcsvtk/isolinepipeline.py index 3406824f0b..4cc1519a03 100644 --- a/Packages/vcs/vcs/vcsvtk/isolinepipeline.py +++ b/Packages/vcs/vcs/vcsvtk/isolinepipeline.py @@ -311,7 +311,8 @@ def _plotInternal(self): z = None kwargs = {"vtk_backend_grid": self._vtkDataSet, "dataset_bounds": self._vtkDataSetBounds, - "plotting_dataset_bounds": plotting_dataset_bounds} + "plotting_dataset_bounds": plotting_dataset_bounds, + "vtk_backend_geo": self._vtkGeoTransform} if ("ratio_autot_viewport" in self._resultDict): kwargs["ratio_autot_viewport"] = vp self._resultDict.update(self._context().renderTemplate( diff --git a/Packages/vcs/vcs/vcsvtk/meshfillpipeline.py b/Packages/vcs/vcs/vcsvtk/meshfillpipeline.py index 64a95c4e31..49320aff93 100644 --- a/Packages/vcs/vcs/vcsvtk/meshfillpipeline.py +++ b/Packages/vcs/vcs/vcsvtk/meshfillpipeline.py @@ -210,7 +210,8 @@ def _plotInternal(self): self._resultDict["vtk_backend_actors"] = actors kwargs = {"vtk_backend_grid": self._vtkDataSet, "dataset_bounds": self._vtkDataSetBounds, - "plotting_dataset_bounds": plotting_dataset_bounds} + "plotting_dataset_bounds": plotting_dataset_bounds, + "vtk_backend_geo": self._vtkGeoTransform} if ("ratio_autot_viewport" in self._resultDict): kwargs["ratio_autot_viewport"] = vp self._template.plot(self._context().canvas, self._data1, self._gm, diff --git a/Packages/vcs/vcs/vcsvtk/vectorpipeline.py b/Packages/vcs/vcs/vcsvtk/vectorpipeline.py index 642884bc6c..c471a6fa93 100644 --- a/Packages/vcs/vcs/vcsvtk/vectorpipeline.py +++ b/Packages/vcs/vcs/vcsvtk/vectorpipeline.py @@ -216,7 +216,8 @@ def _plotInternal(self): create_renderer=True) kwargs = {'vtk_backend_grid': self._vtkDataSet, 'dataset_bounds': self._vtkDataSetBounds, - 'plotting_dataset_bounds': plotting_dataset_bounds} + 'plotting_dataset_bounds': plotting_dataset_bounds, + 'vtk_backend_geo': self._vtkGeoTransform} if ('ratio_autot_viewport' in self._resultDict): kwargs["ratio_autot_viewport"] = vp self._resultDict.update(self._context().renderTemplate( diff --git a/testing/vcs/CMakeLists.txt b/testing/vcs/CMakeLists.txt index bee8b9a45f..b1e6247e72 100644 --- a/testing/vcs/CMakeLists.txt +++ b/testing/vcs/CMakeLists.txt @@ -11,11 +11,26 @@ cdat_add_test(test_vcs_bad_png_path "${PYTHON_EXECUTABLE}" ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_bad_png_path.py ) -cdat_add_test(test_vcs_boxfill_polar - "${PYTHON_EXECUTABLE}" - ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_polar.py - "${BASELINE_DIR}/test_vcs_boxfill_polar.png" - ) + +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} + ) +endforeach() + +foreach(lat_0 45 90) + cdat_add_test(test_vcs_boxfill_orthographic_${lat_0} + "${PYTHON_EXECUTABLE}" + ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_orthographic.py + "${BASELINE_DIR}/test_vcs_boxfill_orthographic_${lat_0}.png" + ${lat_0} + ) +endforeach() + + cdat_add_test(test_vcs_create_get "${PYTHON_EXECUTABLE}" ${cdat_SOURCE_DIR}/testing/vcs/test_vcs_create_get.py diff --git a/testing/vcs/test_vcs_boxfill_orthographic.py b/testing/vcs/test_vcs_boxfill_orthographic.py new file mode 100644 index 0000000000..b0ebbb3a0d --- /dev/null +++ b/testing/vcs/test_vcs_boxfill_orthographic.py @@ -0,0 +1,21 @@ +import os, sys, cdms2, vcs, testing.regression as regression + +baselineName = sys.argv[1] +centerlatitude = float(sys.argv[2]) + + +f = cdms2.open(vcs.sample_data + "/clt.nc") +a = f("clt") + +x = regression.init() +p = x.getprojection('orthographic') +p.centerlatitude = centerlatitude +b = x.createboxfill() +b.projection = p +x.plot(a(latitude=(90,-90)), b, bg=1) + +fileName = os.path.basename(baselineName) +fileName = os.path.splitext(fileName)[0] +fileName += '.png' + +regression.run(x, fileName) diff --git a/testing/vcs/test_vcs_boxfill_polar.py b/testing/vcs/test_vcs_boxfill_projection.py similarity index 64% rename from testing/vcs/test_vcs_boxfill_polar.py rename to testing/vcs/test_vcs_boxfill_projection.py index 869d09802c..6f319efd4b 100644 --- a/testing/vcs/test_vcs_boxfill_polar.py +++ b/testing/vcs/test_vcs_boxfill_projection.py @@ -1,16 +1,20 @@ import os, sys, cdms2, vcs, testing.regression as regression +baselineName = sys.argv[1] +projection = sys.argv[2] + f = cdms2.open(vcs.sample_data + "/clt.nc") a = f("clt") x = regression.init() -p = x.getprojection("polar") +p = x.getprojection(projection) b = x.createboxfill() b.projection = p x.plot(a(latitude=(90,-90)), b, bg=1) -fileName = os.path.basename(__file__) +fileName = os.path.basename(baselineName) fileName = os.path.splitext(fileName)[0] fileName += '.png' -regression.run(x, fileName) \ No newline at end of file + +regression.run(x, fileName)