BUG: proj4 over option causes problems with polar projections

On certain machines, +over option to proj4 results in wrong
projections. It is not clear why there is a difference between different
computers - we have seen this behavior on both unix and mac. +over
does not wrap points outside of -180, 180. To fix those machines we
remove the option for polar projections.

Author: Dan Lipsa

Packages/vcs/Lib/vcs2vtk.py

@@ -570,24 +570,22 @@ def projectArray(w, projection, wc, geo=None):
         pname = projDict.get(projection._type, projection.type)
         projName = pname
         pd.SetName(projName)
-        if projection.type == 'aeqd':
-            # this is a temporary branch to keep the same
-            # baselines
+        if (projection.type == 'aeqd' or
+                projection.type == 'polar stereographic'):
             setProjectionParameters(pd, projection)
+        elif 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:
             pd.SetOptionalParameter("over", "true")
-            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)
-                if ((not hasattr(projection, 'centralmeridian') or
-                     numpy.allclose(projection.centralmeridian, 1e+20))):
-                    pd.SetCentralMeridian(float(xm + xM) / 2.0)
+            setProjectionParameters(pd, projection)
+            if ((not hasattr(projection, 'centralmeridian') or
+                 numpy.allclose(projection.centralmeridian, 1e+20))):
+                pd.SetCentralMeridian(float(xm + xM) / 2.0)
         geo.SetSourceProjection(ps)
         geo.SetDestinationProjection(pd)
 
@@ -613,24 +611,22 @@ def project(pts, projection, wc, geo=None):
         pname = projDict.get(projection._type, projection.type)
         projName = pname
         pd.SetName(projName)
-        if projection.type == 'aeqd':
-            # this is a temporary branch to keep the same
-            # baselines
+        if (projection.type == 'aeqd' or
+                projection.type == 'polar stereographic'):
             setProjectionParameters(pd, projection)
+        elif 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:
             pd.SetOptionalParameter("over", "true")
-            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)
-                if (not hasattr(projection, 'centralmeridian') or
-                        numpy.allclose(projection.centralmeridian, 1e+20)):
-                    pd.SetCentralMeridian(float(xm + xM) / 2.0)
+            setProjectionParameters(pd, projection)
+            if (not hasattr(projection, 'centralmeridian') or
+                    numpy.allclose(projection.centralmeridian, 1e+20)):
+                pd.SetCentralMeridian(float(xm + xM) / 2.0)
         geo.SetSourceProjection(ps)
         geo.SetDestinationProjection(pd)
     geopts = vtk.vtkPoints()