Skip to content

Commit e901bc0

Browse files
danlipsadanlipsa
committed
BUG #1985: orthographic projection plot is empty
This is because proj4 sets points that are not visisble to infinity. We set those points to 0 and hide them.
1 parent 9d06e08 commit e901bc0

File tree

12 files changed

+160
-59
lines changed

12 files changed

+160
-59
lines changed

‎Packages/vcs/vcs/Canvas.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,7 @@
7171
import vcsaddons # noqa
7272
import vcs.manageElements # noqa
7373
import configurator # noqa
74-
from projection import round_projections # noqa
74+
from projection import no_deformation_projections # noqa
7575

7676
# Python < 3 DeprecationWarning ignored by default
7777
warnings.simplefilter('default')
@@ -3497,7 +3497,7 @@ def set_convert_labels(copy_mthd, test=0):
34973497
if hasattr(gm, "priority") and gm.priority == 0:
34983498
return
34993499
p = self.getprojection(gm.projection)
3500-
if p.type in round_projections and (
3500+
if p.type in no_deformation_projections and (
35013501
doratio == "0" or doratio[:4] == "auto"):
35023502
doratio = "1t"
35033503
for keyarg in keyargs.keys():
@@ -3554,7 +3554,7 @@ def set_convert_labels(copy_mthd, test=0):
35543554
t.data.y2 = p.viewport[3]
35553555

35563556
proj = self.getprojection(p.projection)
3557-
if proj.type in round_projections and (
3557+
if proj.type in no_deformation_projections and (
35583558
doratio == "0" or doratio[:4] == "auto"):
35593559
doratio = "1t"
35603560

@@ -3610,7 +3610,7 @@ def set_convert_labels(copy_mthd, test=0):
36103610
tp = "textcombined"
36113611
gm = vcs.elements[tp][arglist[4]]
36123612
p = self.getprojection(gm.projection)
3613-
if p.type in round_projections:
3613+
if p.type in no_deformation_projections:
36143614
doratio = "1t"
36153615
if p.type == 'linear':
36163616
if gm.g_name == 'Gfm':

‎Packages/vcs/vcs/VTKPlots.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -597,6 +597,7 @@ def plot(self, data1, data2, template, gtype, gname, bg, *args, **kargs):
597597

598598
vtk_backend_grid = kargs.get("vtk_backend_grid", None)
599599
vtk_backend_geo = kargs.get("vtk_backend_geo", None)
600+
bounds = vtk_backend_grid.GetBounds() if vtk_backend_grid else None
600601

601602
pipeline = vcsvtk.createPipeline(gm, self)
602603
if pipeline is not None:
@@ -626,7 +627,7 @@ def plot(self, data1, data2, template, gtype, gname, bg, *args, **kargs):
626627
ren,
627628
to=to,
628629
tt=tt,
629-
cmap=self.canvas.colormap)
630+
cmap=self.canvas.colormap, geoBounds=bounds, geo=vtk_backend_geo)
630631
self.setLayer(ren, tt.priority)
631632
self.text_renderers[tt_key] = ren
632633
elif gtype == "line":
@@ -635,7 +636,6 @@ def plot(self, data1, data2, template, gtype, gname, bg, *args, **kargs):
635636
cmap=self.canvas.colormap)
636637
returned["vtk_backend_line_actors"] = actors
637638
create_renderer = True
638-
bounds = vtk_backend_grid.GetBounds() if vtk_backend_grid else None
639639
for act, geo in actors:
640640
ren = self.fitToViewport(
641641
act,

‎Packages/vcs/vcs/projection.py‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -16,13 +16,15 @@
1616
import vcs
1717
import copy
1818

19-
# projection that seems to be doing a circle
20-
# We will probably to add some more in it as we find more that fit this
21-
round_projections = ['polar (non gctp)', 'stereographic',
22-
'orthographic', "ortho", ]
19+
# used to decide if we show longitude labels for round projections or
20+
# latitude labels for elliptical projections
21+
round_projections = ['polar (non gctp)', 'stereographic']
22+
elliptical_projections = ["robinson", "mollweide", 'orthographic', "ortho"]
23+
# projections in this list are not deformed based on the window size
24+
no_deformation_projections = ['polar (non gctp)', 'stereographic',
25+
'orthographic', "ortho", ]
2326

2427
no_over_proj4_parameter_projections = round_projections+["aeqd", "lambert conformal c"]
25-
elliptical_projections = ["robinson", "mollweide"]
2628

2729

2830
def process_src(nm, code):

‎Packages/vcs/vcs/vcs2vtk.py‎

Lines changed: 89 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -4,12 +4,14 @@
44
import numpy
55
import json
66
import os
7+
import math
78
import meshfill
89
from vtk.util import numpy_support as VN
910
import cdms2
1011
import warnings
1112
from projection import round_projections, no_over_proj4_parameter_projections
1213
from vcsvtk import fillareautils
14+
import sys
1315
import numbers
1416

1517
f = open(os.path.join(vcs.prefix, "share", "vcs", "wmo_symbols.json"))
@@ -220,6 +222,34 @@ def getBoundsList(axis, hasCellData, dualGrid):
220222
return None
221223

222224

225+
def setInfToValid(geoPoints, ghost):
226+
'''
227+
Set infinity points to a point that already exists in the list.
228+
If a ghost array is passed, we also hide infinity points.
229+
We return true if any points are infinity
230+
'''
231+
anyInfinity = False
232+
validPoint = [0, 0, 0]
233+
for i in range(geoPoints.GetNumberOfPoints()):
234+
point = geoPoints.GetPoint(i)
235+
if (not math.isinf(point[0]) and not math.isinf(point[1])):
236+
validPoint[0] = point[0]
237+
validPoint[1] = point[1]
238+
break
239+
for i in range(geoPoints.GetNumberOfPoints()):
240+
point = geoPoints.GetPoint(i)
241+
if (math.isinf(point[0]) or math.isinf(point[1])):
242+
anyInfinity = True
243+
newPoint = list(point)
244+
if (math.isinf(point[0])):
245+
newPoint[0] = validPoint[0]
246+
if (math.isinf(point[1])):
247+
newPoint[1] = validPoint[1]
248+
geoPoints.SetPoint(i, newPoint)
249+
ghost.SetValue(i, vtk.vtkDataSetAttributes.HIDDENPOINT)
250+
return anyInfinity
251+
252+
223253
def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False,
224254
dualGrid=False):
225255
continents = False
@@ -444,6 +474,25 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None, genVectors=False,
444474
data1.getAxis(-2))
445475
geo, geopts = project(pts, projection, getWrappedBounds(
446476
wc, [xm, xM, ym, yM], wrap))
477+
# proj4 returns inf for points that are not visible. Set those to a valid point
478+
# and hide them.
479+
ghost = vg.AllocatePointGhostArray()
480+
if (setInfToValid(geopts, ghost)):
481+
# if there are hidden points, we recompute the bounds
482+
xm = ym = sys.float_info.max
483+
xM = yM = - sys.float_info.max
484+
for i in range(pts.GetNumberOfPoints()):
485+
if (ghost.GetValue(i) & vtk.vtkDataSetAttributes.HIDDENPOINT == 0):
486+
# point not hidden
487+
p = pts.GetPoint(i)
488+
if (p[0] < xm):
489+
xm = p[0]
490+
if (p[0] > xM):
491+
xM = p[0]
492+
if (p[1] < ym):
493+
ym = p[1]
494+
if (p[1] > yM):
495+
yM = p[1]
447496
# Sets the vertics into the grid
448497
vg.SetPoints(geopts)
449498
else:
@@ -557,24 +606,42 @@ def apply_proj_parameters(pd, projection, x1, x2, y1, y2):
557606
else:
558607
pd.SetOptionalParameter("over", "false")
559608
setProjectionParameters(pd, projection)
560-
if (hasattr(projection, 'centralmeridian') and
561-
numpy.allclose(projection.centralmeridian, 1e+20)):
562-
pd.SetCentralMeridian(float(x1 + x2) / 2.0)
563-
if (hasattr(projection, 'centerlongitude') and
564-
numpy.allclose(projection.centerlongitude, 1e+20)):
565-
pd.SetOptionalParameter("lon_0", str(float(x1 + x2) / 2.0))
566-
if (hasattr(projection, 'originlatitude') and
567-
numpy.allclose(projection.originlatitude, 1e+20)):
568-
pd.SetOptionalParameter("lat_0", str(float(y1 + y2) / 2.0))
569-
if (hasattr(projection, 'centerlatitude') and
570-
numpy.allclose(projection.centerlatitude, 1e+20)):
571-
pd.SetOptionalParameter("lat_0", str(float(y1 + y2) / 2.0))
572-
if (hasattr(projection, 'standardparallel1') and
573-
numpy.allclose(projection.standardparallel1, 1.e20)):
574-
pd.SetOptionalParameter('lat_1', str(min(y1, y2)))
575-
if (hasattr(projection, 'standardparallel2') and
576-
numpy.allclose(projection.standardparallel2, 1.e20)):
577-
pd.SetOptionalParameter('lat_2', str(max(y1, y2)))
609+
if (hasattr(projection, 'centralmeridian')):
610+
if (numpy.allclose(projection.centralmeridian, 1e+20)):
611+
centralmeridian = float(x1 + x2) / 2.0
612+
else:
613+
centralmeridian = projection.centralmeridian
614+
pd.SetCentralMeridian(centralmeridian)
615+
if (hasattr(projection, 'centerlongitude')):
616+
if (numpy.allclose(projection.centerlongitude, 1e+20)):
617+
centerlongitude = float(x1 + x2) / 2.0
618+
else:
619+
centerlongitude = projection.centerlongitude
620+
pd.SetOptionalParameter("lon_0", str(centerlongitude))
621+
if (hasattr(projection, 'originlatitude')):
622+
if (numpy.allclose(projection.originlatitude, 1e+20)):
623+
originlatitude = float(y1 + y2) / 2.0
624+
else:
625+
originlatitude = projection.originlatitude
626+
pd.SetOptionalParameter("lat_0", str(originlatitude))
627+
if (hasattr(projection, 'centerlatitude')):
628+
if (numpy.allclose(projection.centerlatitude, 1e+20)):
629+
centerlatitude = float(y1 + y2) / 2.0
630+
else:
631+
centerlatitude = projection.centerlatitude
632+
pd.SetOptionalParameter("lat_0", str(centerlatitude))
633+
if (hasattr(projection, 'standardparallel1')):
634+
if (numpy.allclose(projection.standardparallel1, 1.e20)):
635+
standardparallel1 = min(y1, y2)
636+
else:
637+
standardparallel1 = projection.standardparallel1
638+
pd.SetOptionalParameter('lat_1', str(standardparallel1))
639+
if (hasattr(projection, 'standardparallel2')):
640+
if (numpy.allclose(projection.standardparallel2, 1.e20)):
641+
standardparallel2 = max(y1, y2)
642+
else:
643+
standardparallel2 = projection.standardparallel2
644+
pd.SetOptionalParameter('lat_2', str(standardparallel2))
578645

579646

580647
def projectArray(w, projection, wc, geo=None):
@@ -1072,7 +1139,7 @@ def prepTextProperty(p, winSize, to="default", tt="default", cmap=None,
10721139

10731140

10741141
def genTextActor(renderer, string=None, x=None, y=None,
1075-
to='default', tt='default', cmap=None):
1142+
to='default', tt='default', cmap=None, geoBounds=None, geo=None):
10761143
if isinstance(to, str):
10771144
to = vcs.elements["textorientation"][to]
10781145
if isinstance(tt, str):
@@ -1096,21 +1163,8 @@ def genTextActor(renderer, string=None, x=None, y=None,
10961163
sz = renderer.GetRenderWindow().GetSize()
10971164
actors = []
10981165
pts = vtk.vtkPoints()
1099-
geo = None
11001166
if vcs.elements["projection"][tt.projection].type != "linear":
1101-
# Need to figure out new WC
1102-
Npts = 20
1103-
for i in range(Npts + 1):
1104-
X = tt.worldcoordinate[
1105-
0] + float(i) / Npts * (tt.worldcoordinate[1] -
1106-
tt.worldcoordinate[0])
1107-
for j in range(Npts + 1):
1108-
Y = tt.worldcoordinate[
1109-
2] + float(j) / Npts * (tt.worldcoordinate[3] -
1110-
tt.worldcoordinate[2])
1111-
pts.InsertNextPoint(X, Y, 0.)
1112-
geo, pts = project(pts, tt.projection, tt.worldcoordinate, geo=None)
1113-
wc = pts.GetBounds()[:4]
1167+
wc = geoBounds[:4]
11141168
# renderer.SetViewport(tt.viewport[0],tt.viewport[2],tt.viewport[1],tt.viewport[3])
11151169
renderer.SetWorldPoint(wc)
11161170

@@ -1120,8 +1174,8 @@ def genTextActor(renderer, string=None, x=None, y=None,
11201174
prepTextProperty(p, sz, to, tt, cmap)
11211175
pts = vtk.vtkPoints()
11221176
pts.InsertNextPoint(x[i], y[i], 0.)
1123-
if geo is not None:
1124-
geo, pts = project(pts, tt.projection, tt.worldcoordinate, geo=geo)
1177+
if vcs.elements["projection"][tt.projection].type != "linear":
1178+
_, pts = project(pts, tt.projection, tt.worldcoordinate, geo=geo)
11251179
X, Y, tz = pts.GetPoint(0)
11261180
X, Y = world2Renderer(renderer, X, Y, tt.viewport, wc)
11271181
else:

‎Packages/vcs/vcs/vcsvtk/boxfillpipeline.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -152,7 +152,8 @@ def _plotInternal(self):
152152
z = None
153153
kwargs = {"vtk_backend_grid": self._vtkDataSet,
154154
"dataset_bounds": self._vtkDataSetBounds,
155-
"plotting_dataset_bounds": plotting_dataset_bounds}
155+
"plotting_dataset_bounds": plotting_dataset_bounds,
156+
"vtk_backend_geo": self._vtkGeoTransform}
156157
if ("ratio_autot_viewport" in self._resultDict):
157158
kwargs["ratio_autot_viewport"] = vp
158159
self._resultDict.update(self._context().renderTemplate(

‎Packages/vcs/vcs/vcsvtk/isofillpipeline.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -176,7 +176,8 @@ def _plotInternal(self):
176176
z = None
177177
kwargs = {"vtk_backend_grid": self._vtkDataSet,
178178
"dataset_bounds": self._vtkDataSetBounds,
179-
"plotting_dataset_bounds": plotting_dataset_bounds}
179+
"plotting_dataset_bounds": plotting_dataset_bounds,
180+
"vtk_backend_geo": self._vtkGeoTransform}
180181
if ("ratio_autot_viewport" in self._resultDict):
181182
kwargs["ratio_autot_viewport"] = vp
182183
self._resultDict.update(self._context().renderTemplate(

‎Packages/vcs/vcs/vcsvtk/isolinepipeline.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -311,7 +311,8 @@ def _plotInternal(self):
311311
z = None
312312
kwargs = {"vtk_backend_grid": self._vtkDataSet,
313313
"dataset_bounds": self._vtkDataSetBounds,
314-
"plotting_dataset_bounds": plotting_dataset_bounds}
314+
"plotting_dataset_bounds": plotting_dataset_bounds,
315+
"vtk_backend_geo": self._vtkGeoTransform}
315316
if ("ratio_autot_viewport" in self._resultDict):
316317
kwargs["ratio_autot_viewport"] = vp
317318
self._resultDict.update(self._context().renderTemplate(

‎Packages/vcs/vcs/vcsvtk/meshfillpipeline.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -210,7 +210,8 @@ def _plotInternal(self):
210210
self._resultDict["vtk_backend_actors"] = actors
211211
kwargs = {"vtk_backend_grid": self._vtkDataSet,
212212
"dataset_bounds": self._vtkDataSetBounds,
213-
"plotting_dataset_bounds": plotting_dataset_bounds}
213+
"plotting_dataset_bounds": plotting_dataset_bounds,
214+
"vtk_backend_geo": self._vtkGeoTransform}
214215
if ("ratio_autot_viewport" in self._resultDict):
215216
kwargs["ratio_autot_viewport"] = vp
216217
self._template.plot(self._context().canvas, self._data1, self._gm,

‎Packages/vcs/vcs/vcsvtk/vectorpipeline.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -158,7 +158,8 @@ def _plotInternal(self):
158158
create_renderer=True)
159159
kwargs = {'vtk_backend_grid': self._vtkDataSet,
160160
'dataset_bounds': self._vtkDataSetBounds,
161-
'plotting_dataset_bounds': plotting_dataset_bounds}
161+
'plotting_dataset_bounds': plotting_dataset_bounds,
162+
'vtk_backend_geo': self._vtkGeoTransform}
162163
if ('ratio_autot_viewport' in self._resultDict):
163164
kwargs["ratio_autot_viewport"] = vp
164165
self._resultDict.update(self._context().renderTemplate(

‎testing/vcs/CMakeLists.txt‎

Lines changed: 20 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -11,11 +11,26 @@ cdat_add_test(test_vcs_bad_png_path
1111
"${PYTHON_EXECUTABLE}"
1212
${cdat_SOURCE_DIR}/testing/vcs/test_vcs_bad_png_path.py
1313
)
14-
cdat_add_test(test_vcs_boxfill_polar
15-
"${PYTHON_EXECUTABLE}"
16-
${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_polar.py
17-
"${BASELINE_DIR}/test_vcs_boxfill_polar.png"
18-
)
14+
15+
foreach(projection polar mollweide lambert orthographic mercator polyconic robinson)
16+
cdat_add_test(test_vcs_boxfill_${projection}
17+
"${PYTHON_EXECUTABLE}"
18+
${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_projection.py
19+
"${BASELINE_DIR}/test_vcs_boxfill_${projection}.png"
20+
${projection}
21+
)
22+
endforeach()
23+
24+
foreach(lat_0 45 90)
25+
cdat_add_test(test_vcs_boxfill_orthographic_${lat_0}
26+
"${PYTHON_EXECUTABLE}"
27+
${cdat_SOURCE_DIR}/testing/vcs/test_vcs_boxfill_orthographic.py
28+
"${BASELINE_DIR}/test_vcs_boxfill_orthographic_${lat_0}.png"
29+
${lat_0}
30+
)
31+
endforeach()
32+
33+
1934
cdat_add_test(test_vcs_create_get
2035
"${PYTHON_EXECUTABLE}"
2136
${cdat_SOURCE_DIR}/testing/vcs/test_vcs_create_get.py

0 commit comments

Comments
 (0)
⚡