Orthographic

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':

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,

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):

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:

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(

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(

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(

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,

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(

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