possible isofill issue

[doutriaux1]

Celine gave me some data that plot fine in boxfill but seem to be way off in isofill

[aashish24]

@doutriaux1 Can you post the data somewhere?

[doutriaux1]

will do.

[doutriaux1]

@aashish24 I think this is a VTK/Kitware one, I'm doing something wrong but I don't know what.

branch: vcs_celine_isof_bug
test: Test/test_celine_iso_bug.py
On the attached picture, notice how the data around florida is "White" the actual values near there are 1.5 or so, so it should be red-ish
I left the comments on in this branch that shows the lookuptable colors are set properly as well as the levels. Not sure why it is white. "sFilter" issue? Bandedcontour issue?

[doutriaux1]

Also on this picture you can see why jerry thinks it looks "horrible" everythig is indeed "blurry" (look at the red colors). bug #512

[aashish24]

@doutriaux1 are you converting the data to the image somehow? The polydata cannot be blurry. That would be really odd. I guess we have to look at your code together?

[doutriaux1]

using a BUNCH of filter
sFilter = vtk.vtkDataSetSurfaceFilter()
sFilter.SetInputData(ug) # ug here is a structured grid
mapper = vtk.vtkPolyDataMapper()
cot = vtk.vtkBandedPolyDataContourFilter()
cot.ClippingOn()
cot.SetInputData(sFilter.GetOutput())
cot.SetNumberOfContours(len(l))
for j,v in enumerate(l):
cot.SetValue(j,l[j])
cot.Update()
lut = vtk.vtkLookupTable()
print "NCOLORS:",len(COLS[i])
lut.SetNumberOfTableValues(len(COLS[i]))
for j,color in enumerate(COLS[i]):
r,g,b = cmap.index[color]
print l[j],vcs.colors.rgb2str(r_2.55,g_2.55,b*2.55),l[j+1]
lut.SetTableValue(j,r/100.,g/100.,b/100.)
mapper.SetLookupTable(lut)
mapper.SetScalarRange(lmn,lmx)

[doutriaux1]

will post here a pure vtk code that does what I do. it wil be easier to debug

[aashish24]

@doutriaux1 thanks.

[doutriaux1]

here you go:

from vtk.util import numpy_support as VN
f=cdms2.open("celine.nc")
data=f("data")
ug = vtk.vtkStructuredGrid()
lat = data.getLatitude()
lon = data.getLongitude()
lat = lat[:,numpy.newaxis]_numpy.ones(lon.shape)[numpy.newaxis,:]
lon = lon[numpy.newaxis,:]_numpy.ones(lat.shape)[:,numpy.newaxis]
ug.SetDimensions(lat.shape[1],lat.shape[0],1)
lon = numpy.ma.ravel(lon)
lat = numpy.ma.ravel(lat)
sh = list(lat.shape)
sh.append(1)
lon = numpy.ma.reshape(lon,sh)
lat = numpy.ma.reshape(lat,sh)
z = numpy.zeros(lon.shape)
m3 = numpy.concatenate((lon,lat),axis=1)
m3 = numpy.concatenate((m3,z),axis=1)
xm=lon.min()
xM=lon.max()
ym=lat.min()
yM=lat.max()

# First create the points/vertices (in vcs terms)

deep = True
pts = vtk.vtkPoints()

## Convert nupmy array to vtk ones

ppV = VN.numpy_to_vtk(m3,deep=deep)
pts.SetData(ppV)
ug.SetPoints(pts)
data = VN.numpy_to_vtk(data.filled(0.).flat,deep=True)
ug.GetPointData().SetScalars(data)
sFilter = vtk.vtkDataSetSurfaceFilter()
sFilter.SetInputData(ug)
sFilter.Update()
l = [-4.0584776401519775, -2, -1.5, -1.0, -0.4, -0.1, 0, 0.1, 0.4, 1.0, 1.5, 2, 7.875967502593994]
cols = [[0, 57, 100],[0, 100, 33],[27, 100, 0],[88, 100, 0], [100, 78, 0], [100, 100, 100], [100, 100, 100], [100, 32, 0], [100, 17, 0], [75, 0, 0], [45, 0, 0], [49, 0, 24],]
mapper = vtk.vtkPolyDataMapper()
cot = vtk.vtkBandedPolyDataContourFilter()
cot.ClippingOn()
cot.SetInputData(sFilter.GetOutput())
cot.SetNumberOfContours(len(l))
for j,v in enumerate(l):
  cot.SetValue(j,l[j])
cot.Update()
mapper.SetInputConnection(cot.GetOutputPort())
lut = vtk.vtkLookupTable()
lut.SetNumberOfTableValues(len(cols))
for j,color in enumerate(cols):
  r,g,b = cols[j]
  lut.SetTableValue(j,r/100.,g/100.,b/100.)
mapper.SetLookupTable(lut)
mapper.SetScalarRange(-4.05847764015,7.87596750259)
act = vtk.vtkActor()
act.SetMapper(mapper)
ren=vtk.vtkRenderer()
ren.AddActor(act)
renWin=vtk.vtkRenderWindow()
renWin.SetSize(800,600)
renWin.AddRenderer(ren)
renWin.Render()
raw_input("ok press enter now")

[aashish24]

Ok. @jbeezley and I are looking to the code now. We will report back once we know what's going on.

[doutriaux1]

Ok thank you, thank you . thank you! I suspect (but no idea) a celltodata or datatocell call needs (not) to be made somewhere. BTW Plotting in boxfill the data comes out correctly

[jbeezley]

I'm far from a VTK expert, but I made a few changes that seem to help.

import vtk,cdms2,numpy

from vtk.util import numpy_support as VN

f=cdms2.open("celine.nc")
data=f("data")
ug = vtk.vtkStructuredGrid()
lat = data.getLatitude()
lon = data.getLongitude()
lat = lat[:,numpy.newaxis]*numpy.ones(lon.shape)[numpy.newaxis,:]
lon = lon[numpy.newaxis,:]*numpy.ones(lat.shape)[:,numpy.newaxis]
ug.SetDimensions(lat.shape[1],lat.shape[0],1)
lon = numpy.ma.ravel(lon)
lat = numpy.ma.ravel(lat)
sh = list(lat.shape)
sh.append(1)
lon = numpy.ma.reshape(lon,sh)
lat = numpy.ma.reshape(lat,sh)
z = numpy.zeros(lon.shape)
m3 = numpy.concatenate((lon,lat),axis=1)
m3 = numpy.concatenate((m3,z),axis=1)
xm=lon.min()
xM=lon.max()
ym=lat.min()
yM=lat.max()
# First create the points/vertices (in vcs terms)
deep = True
pts = vtk.vtkPoints()
## Convert nupmy array to vtk ones
ppV = VN.numpy_to_vtk(m3,deep=deep)
pts.SetData(ppV)
ug.SetPoints(pts)
data = VN.numpy_to_vtk(data.filled(0.).flat,deep=True)
ug.GetPointData().SetScalars(data)

sFilter = vtk.vtkDataSetSurfaceFilter()
sFilter.SetInputData(ug)
sFilter.Update()
l = [-4.0584776401519775, -2, -1.5, -1.0, -0.4, -0.1, 0, 0.1, 0.4, 1.0, 1.5, 2, 7.875967502593994]
cols = [[0, 57, 100],[0, 100, 33],[27, 100, 0],[88, 100, 0], [100, 78, 0], [100, 100, 100], [100, 100, 100], [100, 32, 0], [100, 17, 0], [75, 0, 0], [45, 0, 0], [49, 0, 24],]

mapper = vtk.vtkPolyDataMapper()
cot = vtk.vtkBandedPolyDataContourFilter()
cot.ClippingOn()
cot.SetInputData(sFilter.GetOutput())
cot.SetNumberOfContours(len(l))
for j,v in enumerate(l):
  cot.SetValue(j,l[j])
cot.SetScalarModeToIndex()
cot.Update()

mapper.SetInputConnection(cot.GetOutputPort())
lut = vtk.vtkLookupTable()
lut.SetNumberOfTableValues(len(cols))
for j,color in enumerate(cols):
  r,g,b = cols[j]
  lut.SetTableValue(j,r/100.,g/100.,b/100.)
mapper.SetLookupTable(lut)
mapper.SetScalarRange(0, len(l) - 1)
mapper.SetScalarModeToUseCellData()

act = vtk.vtkActor()
act.SetMapper(mapper)

ren=vtk.vtkRenderer()
ren.AddActor(act)
renWin=vtk.vtkRenderWindow()
renWin.SetSize(800,600)
renWin.AddRenderer(ren)
renWin.Render()

iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
ren.SetBackground(1, 1, 1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
iren.Start()
raw_input("ok press enter now")

The contours need to be adjusted, but here is the result.

[doutriaux1]

yes! I confirmed it works. AS I suspected it had something to do with celldata vs pointdata.

@aashish24 since we set the data on point, why do we need to use:
mapper.SetScalarModeToUseCellData()