Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BUG #1811: Show point information for plots using geographic projecti… #1878

Merged
merged 2 commits into from Mar 21, 2016
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
166 changes: 112 additions & 54 deletions Packages/vcs/Lib/VTKPlots.py
Expand Up @@ -45,18 +45,27 @@ def __init__(self, canvas, renWin=None, debug=False, bg=None, geometry=None):
self.renderer = None
self._renderers = {}
self._plot_keywords = [
'renderer',
'vtk_backend_grid',
'cdmsfile',
'cell_coordinates'
# used to render the continents
'continents_renderer',
# dataset bounds in lon/lat coordinates
'dataset_bounds',
# the same as vcs.utils.getworldcoordinates for now. getworldcoordinates uses
# gm.datawc_... or, if that is not set, it uses data axis margins (without bounds).
'plotting_dataset_bounds',
# this may be slightly smaller than the data viewport. It is used
# for vectorpipeline when drawn on top of a boxfill for instance
'dataset_viewport',
# used to render the dataset
'dataset_renderer',
# dataset scale: (xScale, yScale)
'dataset_scale',
# the same as vcs.utils.getworldcoordinates for now. getworldcoordinates uses
# gm.datawc_... or, if that is not set, it uses data axis margins (without bounds).
'plotting_dataset_bounds',
'renderer',
'vtk_backend_grid',
# vtkGeoTransform used for geographic transformation
'vtk_backend_geo',
'cdmsfile',
'cell_coordinates']
]
self.numberOfPlotCalls = 0
self.renderWindowSize = None
self.clickRenderer = None
Expand Down Expand Up @@ -131,47 +140,97 @@ def leftButtonPressEvent(self, obj, event):
continue
t = vcs.elements["template"][d.template]
gm = vcs.elements[d.g_type][d.g_name]
if gm.projection != "linear":
return
if t.data.x1 <= x <= t.data.x2 and t.data.y1 <= y <= t.data.y2:
x1, x2, y1, y2 = vcs.utils.getworldcoordinates(gm,
d.array[0].getAxis(-1),
d.array[0].getAxis(-2))

X = (x - t.data.x1) / (t.data.x2 - t.data.x1) * (x2 - x1) + x1
Y = (y - t.data.y1) / (t.data.y2 - t.data.y1) * (y2 - y1) + y1

# Ok we now have the X/Y values we need to figure out the
# indices
try:
I = d.array[0].getAxis(-1).mapInterval((X, X, 'cob'))[0]
# for non-linear projection or for meshfill. Meshfill is wrapped at
# VTK level, so vcs calculations do not work.
if gm.projection != "linear" or gm.g_name == 'Gfm':
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@danlipsa why not use the same code (under if) for both cases?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because this code does not work now for isofill, isoline and vectors. When we fix dataset creation we can get rid of the vcs code.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@danlipsa thanks. And I guess we have an issue for it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes we do:
#1881

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks @danlipsa

selector = vtk.vtkHardwareSelector()
datasetRenderer = d.backend['dataset_renderer']
continentsRenderer = d.backend.get('continents_renderer')
dataset = d.backend['vtk_backend_grid']
if (datasetRenderer and dataset):
selector.SetRenderer(datasetRenderer)
selector.SetArea(xy[0], xy[1], xy[0], xy[1])
selector.SetFieldAssociation(vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS)
# We want to be able see information behind continents
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@danlipsa why field association?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because we want cellids which will give us cell scalars.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I see thanks. So are we not looking for point scalars in this case?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is correct. For point data, the user has to click fairly close the a point to see the value. I have to experiment with that, I am not sure.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is correct. For point data, the user has to click fairly close the a point to see the value. I have to experiment with that, I am not sure

Do we have test for it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not implemented yet. I'll do it together with or after
#1881

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not implemented yet. I'll do it together with or after
#1881

thanks. Should we file an issue for it @danlipsa ?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if (continentsRenderer):
continentsRenderer.SetDraw(False)
selection = selector.Select()
if (continentsRenderer):
continentsRenderer.SetDraw(True)
if (selection.GetNumberOfNodes() > 0):
selectionNode = selection.GetNode(0)
prop = selectionNode.GetProperties().Get(selectionNode.PROP())
if (prop):
cellIds = prop.GetMapper().GetInput().GetCellData().GetGlobalIds()
if (cellIds):
# scalar value
a = selectionNode.GetSelectionData().GetArray(0)
geometryId = a.GetValue(0)
cellId = cellIds.GetValue(geometryId)
scalars = dataset.GetCellData().GetScalars()
value = scalars.GetValue(cellId)
geoTransform = d.backend['vtk_backend_geo']
if (geoTransform):
geoTransform.Inverse()
# Use the world picker to get world coordinates
# we deform the dataset, so we need to fix the
# world picker using xScale, yScale
xScale, yScale = d.backend['dataset_scale']
worldPicker = vtk.vtkWorldPointPicker()
worldPicker.Pick(xy[0], xy[1], 0, datasetRenderer)
worldPosition = list(worldPicker.GetPickPosition())
if (xScale > yScale):
worldPosition[0] /= (xScale/yScale)
else:
worldPosition[1] /= (yScale/xScale)
lonLat = worldPosition
if (geoTransform):
geoTransform.InternalTransformPoint(worldPosition, lonLat)
geoTransform.Inverse()
st += "Var: %s\n" % d.array[0].id
if (float("inf") not in lonLat):
st += "X=%4.1f\nY=%4.1f\n" % (lonLat[0], lonLat[1])
st += "Value: %g" % value
else:
if t.data.x1 <= x <= t.data.x2 and t.data.y1 <= y <= t.data.y2:
x1, x2, y1, y2 = vcs.utils.getworldcoordinates(gm,
d.array[0].getAxis(-1),
d.array[0].getAxis(-2))

X = (x - t.data.x1) / (t.data.x2 - t.data.x1) * (x2 - x1) + x1
Y = (y - t.data.y1) / (t.data.y2 - t.data.y1) * (y2 - y1) + y1

# Ok we now have the X/Y values we need to figure out the
# indices
try:
J = d.array[
0].getAxis(-2).mapInterval((Y, Y, 'cob'))[0]
# Values at that point
V = d.array[0][..., J, I]
except:
V = d.array[0][..., I]
if isinstance(V, numpy.ndarray):
# Grab the appropriate time slice
if self.canvas.animate.created():
t = self.canvas.animate.frame_num
try:
taxis = V.getTime()
V = V(time=taxis[t % len(taxis)]).flat[0]
except:
I = d.array[0].getAxis(-1).mapInterval((X, X, 'cob'))[0]
try:
J = d.array[
0].getAxis(-2).mapInterval((Y, Y, 'cob'))[0]
# Values at that point
V = d.array[0][..., J, I]
except:
V = d.array[0][..., I]
if isinstance(V, numpy.ndarray):
# Grab the appropriate time slice
if self.canvas.animate.created():
t = self.canvas.animate.frame_num
try:
taxis = V.getTime()
V = V(time=taxis[t % len(taxis)]).flat[0]
except:
V = V.flat[0]
else:
V = V.flat[0]
else:
V = V.flat[0]
try:
st += "Var: %s\nX[%i] = %g\nY[%i] = %g\nValue: %g" % (
d.array[0].id, I, X, J, Y, V)
try:
st += "Var: %s\nX[%i] = %4.1f\nY[%i] = %4.1f\nValue: %g" % (
d.array[0].id, I, X, J, Y, V)
except:
st += "Var: %s\nX = %4.1f\nY[%i] = %4.1f\nValue: %g" % (
d.array[0].id, X, I, Y, V)
except:
st += "Var: %s\nX = %g\nY[%i] = %g\nValue: %g" % (
d.array[0].id, X, I, Y, V)
except:
st += "Var: %s\nX=%g\nY=%g\nValue = N/A" % (
d.array[0].id, X, Y)
st += "Var: %s\nX=%g\nY=%g\nValue = N/A" % (
d.array[0].id, X, Y)
if st == "":
return
ren = vtk.vtkRenderer()
Expand Down Expand Up @@ -740,13 +799,12 @@ def plotContinents(self, wc, projection, wrap, vp, priority, **kargs):
# Stippling
vcs2vtk.stippleLine(line_prop, contLine.type[0])
vtk_backend_grid = kargs.get("vtk_backend_grid", None)
self.fitToViewportBounds(contActor,
vp,
wc=wc, geo=geo,
geoBounds=vtk_backend_grid.GetBounds(),
priority=priority,
create_renderer=True)
return {}
return self.fitToViewportBounds(contActor,
vp,
wc=wc, geo=geo,
geoBounds=vtk_backend_grid.GetBounds(),
priority=priority,
create_renderer=True)

def renderTemplate(self, tmpl, data, gm, taxis, zaxis, **kargs):
# ok first basic template stuff, let's store the displays
Expand Down Expand Up @@ -1199,7 +1257,7 @@ def fitToViewportBounds(self, Actor, vp, wc=None, geoBounds=None, geo=None, prio

# Data range in World Coordinates
if priority == 0:
return None
return (None, 1, 1)
vp = tuple(vp)
if wc is None:
Xrg = list(Actor.GetXRange())
Expand Down Expand Up @@ -1315,7 +1373,7 @@ def fitToViewportBounds(self, Actor, vp, wc=None, geoBounds=None, geo=None, prio
plane = planeCollection.GetNextItem()

Renderer.AddActor(Actor)
return Renderer
return (Renderer, xScale, yScale)

def update_input(self, vtkobjects, array1, array2=None, update=True):
if "vtk_backend_grid" in vtkobjects:
Expand Down
13 changes: 12 additions & 1 deletion Packages/vcs/Lib/vcs2vtk.py
Expand Up @@ -507,6 +507,10 @@ def genGrid(data1, data2, gm, deep=True, grid=None, geo=None):
vg.SetPoints(geopts)
else:
vg = grid
# Add a GlobalIds array to keep track of cell ids throughout the pipeline
globalIds = numpy_to_vtk_wrapper(numpy.arange(0, vg.GetNumberOfCells()), deep=True)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@danlipsa why we have to do this? Meaning why we have to explicitly append the global ids?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because the hardware selector renders and gets the cell ids of the geometry which are not usually the same as the cell ids of the dataset.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because the hardware selector renders and gets the cell ids of the geometry which are not usually the same as the cell ids of the dataset.

this is interesting. But by default should they be same?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do you mean 'by default'? In the uvcdat case we have threshold + different mapper for different scalar ranges so the ids are not the dataset ids. If you have only the dataset in the pipeline I think only the polydata will not change ids. Maybe there are others but I would not bet on it :-)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you have only the dataset in the pipeline I think only the polydata will not change ids. Maybe there are others but I would not bet on it :-)

I see. yeah, I won't either.

globalIds.SetName('GlobalIds')
vg.GetCellData().SetGlobalIds(globalIds)
out = {"vtk_backend_grid": vg,
"xm": xm,
"xM": xM,
Expand Down Expand Up @@ -867,6 +871,9 @@ def doWrap(Act, wc, wrap=[0., 360], fastClip=True):
return Act
Mapper = Act.GetMapper()
data = Mapper.GetInput()
# insure that GLOBALIDS are not removed by the append filter
attributes = data.GetCellData()
attributes.SetActiveAttribute(-1, attributes.GLOBALIDS)
xmn = min(wc[0], wc[1])
xmx = max(wc[0], wc[1])
if numpy.allclose(xmn, 1.e20) or numpy.allclose(xmx, 1.e20):
Expand All @@ -878,7 +885,6 @@ def doWrap(Act, wc, wrap=[0., 360], fastClip=True):
ymx = abs(wrap[0])
ymn = -wrap[0]

# Prepare MultiBlock and puts in oriinal data
appendFilter = vtk.vtkAppendPolyData()
appendFilter.AddInputData(data)
appendFilter.Update()
Expand Down Expand Up @@ -954,6 +960,11 @@ def doWrap(Act, wc, wrap=[0., 360], fastClip=True):
clipper.SetClipFunction(clipBox)
clipper.SetInputConnection(appendFilter.GetOutputPort())
clipper.Update()
# set globalids attribute
attributes = clipper.GetOutput().GetCellData()
globalIdsIndex = vtk.mutable(-1)
attributes.GetArray("GlobalIds", globalIdsIndex)
attributes.SetActiveAttribute(globalIdsIndex, attributes.GLOBALIDS)

Mapper.SetInputData(clipper.GetOutput())
return Act
Expand Down
27 changes: 14 additions & 13 deletions Packages/vcs/Lib/vcsvtk/boxfillpipeline.py
Expand Up @@ -96,12 +96,7 @@ def _updateContourLevelsAndColorsForBoxfill(self):
def _createPolyDataFilter(self):
"""Overrides baseclass implementation."""
self._vtkPolyDataFilter = vtk.vtkDataSetSurfaceFilter()
if self._useCellScalars:
p2c = vtk.vtkPointDataToCellData()
p2c.SetInputData(self._vtkDataSet)
self._vtkPolyDataFilter.SetInputConnection(p2c.GetOutputPort())
else:
self._vtkPolyDataFilter.SetInputData(self._vtkDataSet)
self._vtkPolyDataFilter.SetInputData(self._vtkDataSet)

def _plotInternal(self):
"""Overrides baseclass implementation."""
Expand All @@ -126,6 +121,8 @@ def _plotInternal(self):
_style = self._gm.fillareastyle
vp = [self._template.data.x1, self._template.data.x2,
self._template.data.y1, self._template.data.y2]
dataset_renderer = None
xScale, yScale = (1, 1)
for mapper in self._mappers:
act = vtk.vtkActor()
act.SetMapper(mapper)
Expand Down Expand Up @@ -170,12 +167,14 @@ def _plotInternal(self):

# create a new renderer for this mapper
# (we need one for each mapper because of camera flips)
self._context().fitToViewportBounds(
dataset_renderer, xScale, yScale = self._context().fitToViewportBounds(
act, vp,
wc=plotting_dataset_bounds, geoBounds=self._vtkDataSet.GetBounds(),
geo=self._vtkGeoTransform,
priority=self._template.data.priority,
create_renderer=True)
create_renderer=(dataset_renderer is None))
self._resultDict['dataset_renderer'] = dataset_renderer
self._resultDict['dataset_scale'] = (xScale, yScale)

for act in patternActors:
if self._vtkGeoTransform is None:
Expand Down Expand Up @@ -247,11 +246,13 @@ def _plotInternal(self):
self._useContinents = False
if self._useContinents:
projection = vcs.elements["projection"][self._gm.projection]
self._context().plotContinents(plotting_dataset_bounds, projection,
self._dataWrapModulo,
vp, self._template.data.priority,
vtk_backend_grid=self._vtkDataSet,
dataset_bounds=self._vtkDataSetBounds)
continents_renderer, xScale, yScale = self._context().plotContinents(
plotting_dataset_bounds, projection,
self._dataWrapModulo,
vp, self._template.data.priority,
vtk_backend_grid=self._vtkDataSet,
dataset_bounds=self._vtkDataSetBounds)
self._resultDict['continents_renderer'] = continents_renderer

def _plotInternalBoxfill(self):
"""Implements the logic to render a non-custom boxfill."""
Expand Down
21 changes: 13 additions & 8 deletions Packages/vcs/Lib/vcsvtk/isofillpipeline.py
Expand Up @@ -136,6 +136,8 @@ def _plotInternal(self):
ct = 0
vp = [self._template.data.x1, self._template.data.x2,
self._template.data.y1, self._template.data.y2]
dataset_renderer = None
xScale, yScale = (1, 1)
for mapper in mappers:
act = vtk.vtkActor()
act.SetMapper(mapper)
Expand Down Expand Up @@ -174,13 +176,14 @@ def _plotInternal(self):

# create a new renderer for this mapper
# (we need one for each mapper because of cmaera flips)
self._context().fitToViewportBounds(
dataset_renderer, xScale, yScale = self._context().fitToViewportBounds(
act, vp,
wc=plotting_dataset_bounds, geoBounds=self._vtkDataSet.GetBounds(),
geo=self._vtkGeoTransform,
priority=self._template.data.priority,
create_renderer=True)

create_renderer=(dataset_renderer is None))
self._resultDict['dataset_renderer'] = dataset_renderer
self._resultDict['dataset_scale'] = (xScale, yScale)
for act in patternActors:
self._context().fitToViewportBounds(
act, vp,
Expand Down Expand Up @@ -241,8 +244,10 @@ def _plotInternal(self):
self._useContinents = False
if self._useContinents:
projection = vcs.elements["projection"][self._gm.projection]
self._context().plotContinents(plotting_dataset_bounds, projection,
self._dataWrapModulo,
vp, self._template.data.priority,
vtk_backend_grid=self._vtkDataSet,
dataset_bounds=self._vtkDataSetBounds)
continents_renderer, xScale, yScale = self._context().plotContinents(
plotting_dataset_bounds, projection,
self._dataWrapModulo,
vp, self._template.data.priority,
vtk_backend_grid=self._vtkDataSet,
dataset_bounds=self._vtkDataSetBounds)
self._resultDict['continents_renderer'] = continents_renderer
21 changes: 13 additions & 8 deletions Packages/vcs/Lib/vcsvtk/isolinepipeline.py
Expand Up @@ -156,6 +156,8 @@ def _plotInternal(self):
countLevels = 0
vp = [self._template.data.x1, self._template.data.x2,
self._template.data.y1, self._template.data.y2]
dataset_renderer = None
xScale, yScale = (1, 1)
for i, l in enumerate(tmpLevels):
numLevels = len(l)

Expand Down Expand Up @@ -289,15 +291,16 @@ def _plotInternal(self):

# create a new renderer for this mapper
# (we need one for each mapper because of cmaera flips)
self._context().fitToViewportBounds(
dataset_renderer, xScale, yScale = self._context().fitToViewportBounds(
act, vp,
wc=plotting_dataset_bounds, geoBounds=self._vtkDataSet.GetBounds(),
geo=self._vtkGeoTransform,
priority=self._template.data.priority,
create_renderer=True)
create_renderer=(dataset_renderer is None))

countLevels += len(l)

self._resultDict['dataset_renderer'] = dataset_renderer
self._resultDict['dataset_scale'] = (xScale, yScale)
if len(textprops) > 0:
self._resultDict["vtk_backend_contours_labels_text_properties"] = \
textprops
Expand Down Expand Up @@ -348,8 +351,10 @@ def _plotInternal(self):
self._useContinents = False
if self._useContinents:
projection = vcs.elements["projection"][self._gm.projection]
self._context().plotContinents(plotting_dataset_bounds, projection,
self._dataWrapModulo,
vp, self._template.data.priority,
vtk_backend_grid=self._vtkDataSet,
dataset_bounds=self._vtkDataSetBounds)
continents_renderer, xScale, yScale = self._context().plotContinents(
plotting_dataset_bounds, projection,
self._dataWrapModulo,
vp, self._template.data.priority,
vtk_backend_grid=self._vtkDataSet,
dataset_bounds=self._vtkDataSetBounds)
self._resultDict['continents_renderer'] = continents_renderer