Main Page | Report Page

 

  Computers Forum Index » Computer Languages (IDL-PVWAVE) » problem with TRIANGULATION option in CONTOUR...

Author Message
Ardhuin...
Posted: Sat Oct 23, 2010 12:35 pm
 
Dear all,
I have been having problems with plotting output of a numerical model
that uses unstructured grids using IDL: this model computes wave
heights in the ocean. If I use the set of triangles from the model, I
get really funny errors:
Out of range subscript encountered: <LONG Array[66453]>.
although I only have 9626 points and about 16000 triangles.

The command I use is
CONTOUR,tablep,x,y,$
xstyle=5,ystyle=5,/FOLLOW,/CELL_FILL, TRIANGULATION=tri,
$
C_COLOR=colorind(0:c_numlev-2+addmini+addmaxi), $
LEVELS=lev,/
NOERASE,TITLE=title,CLIP=[rangex(0),rangey(0),rangex(1),rangey(1)], $
XRANGE=rangex,YRANGE=rangey,MAX_VALUE=maxval, $
POSITION=[blx*winx/mwinx,bly*winy/mwiny,trx*winx/
mwinx,try*winy/mwiny]

Where tablep , x and y
If I first do TRIANGULATE,X,Y,tri then the contours comes out OK...
but they cut out through land boundaries and islands.

So I was thinking: my triangles must be wrong ...
but if I do a TRIGRID with my triangles then I can plot nicely with
TV ...
array=trigrid(X,Y,tablep,tri,[dx,dy],
[rangex(0),rangey(0),rangex(1),rangey(1)], $
MAX_VALUE=maxval)

So my triangles are OK for TRIGRID but not for CONTOUR... How is that
possible ??
 
David Fanning...
Posted: Sat Oct 23, 2010 4:57 pm
 
Ardhuin writes:

Quote:
So my triangles are OK for TRIGRID but not for CONTOUR... How is that
possible ??

Have you tried your contour plot without the CLIP keyword?
What is CLIP doing for you?

Cheers,

David


--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
 
Ardhuin...
Posted: Sat Oct 23, 2010 10:00 pm
 
Yes, the CLIP was useless, but even simplifying the command to

CONTOUR,tablep,x,y,$
xstyle=5,ystyle=5,/FOLLOW,/CELL_FILL, TRIANGULATION=tri, $
C_COLOR=colorind(0:c_numlev-2+addmini+addmaxi), $
LEVELS=lev,/
XRANGE=rangex,YRANGE=rangey, $
POSITION=[blx*winx/mwinx,bly*winy/mwiny,trx*winx/
mwinx,try*winy/mwiny]

I do get the same error...

Fabrice Ardhuin
 
Ben Tupper...
Posted: Sun Oct 24, 2010 2:52 am
 
On 10/23/10 8:35 AM, Ardhuin wrote:
Quote:
Dear all,
I have been having problems with plotting output of a numerical model
that uses unstructured grids using IDL: this model computes wave
heights in the ocean. If I use the set of triangles from the model, I
get really funny errors:
Out of range subscript encountered:<LONG Array[66453]>.
although I only have 9626 points and about 16000 triangles.

The command I use is
CONTOUR,tablep,x,y,$
xstyle=5,ystyle=5,/FOLLOW,/CELL_FILL, TRIANGULATION=tri,
$
C_COLOR=colorind(0:c_numlev-2+addmini+addmaxi), $
LEVELS=lev,/
NOERASE,TITLE=title,CLIP=[rangex(0),rangey(0),rangex(1),rangey(1)], $
XRANGE=rangex,YRANGE=rangey,MAX_VALUE=maxval, $
POSITION=[blx*winx/mwinx,bly*winy/mwiny,trx*winx/
mwinx,try*winy/mwiny]

Where tablep , x and y
If I first do TRIANGULATE,X,Y,tri then the contours comes out OK...
but they cut out through land boundaries and islands.

So I was thinking: my triangles must be wrong ...
but if I do a TRIGRID with my triangles then I can plot nicely with
TV ...
array=trigrid(X,Y,tablep,tri,[dx,dy],
[rangex(0),rangey(0),rangex(1),rangey(1)], $
MAX_VALUE=maxval)

So my triangles are OK for TRIGRID but not for CONTOUR... How is that
possible ??


Hi,

You might consider running the points through GRID_INPUT first. I
really have no idea if that will help in this case, but it often comes
to the rescue when interpolating.

Cheers,
Ben
 
David Fanning...
Posted: Sun Oct 24, 2010 3:15 am
 
Ben Tupper writes:

Quote:
You might consider running the points through GRID_INPUT first. I
really have no idea if that will help in this case, but it often comes
to the rescue when interpolating.

Yes, that's a good suggestion.

Or, I would try just gridding the data with Triangulate
and Trigrid, and pass the gridded data to Contour.
In other words, don't try to use the TRIANGULATION
keyword. (I have to admit, I had never heard of this
keyword!)

Here are the commands I use to grid irregular data
from a book that will be hot of the press in a month
or so. I get the X and Y vectors that go with the
gridded data from the XGRID and YGRID output
keywords to TriGrid.

Triangulate, lonIrr, latIrr, triangles
gridData = Trigrid(lonIrr, latIrr, dataIrr, $
triangles, NX=41, NY=41, $
XGrid=xgrid, YGrid=ygrid)
Contour, gridData, xgrid, ygrid, /Cell_Fill, $
Levels=levels, Background=FSC_Color('white'),$
Position=[0.125, 0.125, 0.95, 0.80], $
Color=FSC_Color('black'), XStyle=1, YStyle=1, $
C_Colors=Indgen(nlevels)+1

Cheers,

David

P.S. If you want to send the data, I'm always looking
for perverse examples to add to the book! :-)

--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
 
Ardhuin...
Posted: Mon Oct 25, 2010 7:06 pm
 
Well, this is exactly what I do not want to do: I want to use my
triangulation because it contains the information on the grid
connectivity (islands are "holes" in the grid where the contouring
should not do any filling).

I have put the basic dataset at this ftp address:
ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/wavewatch3/TOOLS/IDL/BUG/
with some illustrative plots. The standard contour fills everything,
including land. If I do a TRIGRID before a TV then I get want I want:
land and Islands are blank... but I'd like to do this with contour to
use the native data resolution...

Fabrice
 
David Fanning...
Posted: Tue Oct 26, 2010 2:17 am
 
Ardhuin writes:

Quote:
Well, this is exactly what I do not want to do: I want to use my
triangulation because it contains the information on the grid
connectivity (islands are "holes" in the grid where the contouring
should not do any filling).

I have put the basic dataset at this ftp address:
ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/wavewatch3/TOOLS/IDL/BUG/
with some illustrative plots. The standard contour fills everything,
including land. If I do a TRIGRID before a TV then I get want I want:
land and Islands are blank... but I'd like to do this with contour to
use the native data resolution...

It seemed to me you wanted to put a contour plot on top
of an image, so that's what I've shown you how to do.
You will need Coyote Library programs to run the code.

http://www.dfanning.com/programs/coyoteprograms.zip

You can find the program I wrote, and a PNG file of what
the output looks like here:

http://www.dfanning.com/misc/contour_on_image.zip

Cheers,

David


--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
 
Ben Tupper...
Posted: Tue Oct 26, 2010 4:29 am
 
On 10/25/10 3:06 PM, Ardhuin wrote:
Quote:
Well, this is exactly what I do not want to do: I want to use my
triangulation because it contains the information on the grid
connectivity (islands are "holes" in the grid where the contouring
should not do any filling).

I have put the basic dataset at this ftp address:
ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/wavewatch3/TOOLS/IDL/BUG/
with some illustrative plots. The standard contour fills everything,
including land. If I do a TRIGRID before a TV then I get want I want:
land and Islands are blank... but I'd like to do this with contour to
use the native data resolution...


Hi,

When I load the data you have posted and generate the triangulation I
get a different number of triangles than you. You can see this in the
output of help - note that tri2 has 19160 vertices. Also, you clip the
maximum Z value when you use TRIGRID but you don't specify it for
CONTOUR. I don't know if these are important issues, but perhaps you
are unwittingly comparing apples to oranges?

Cheers,
Ben

IDL> .full_reset
IDL> restore, filename =
"/Users/Shared/data/tri/TRIANGULATION_MARSEILLE.sav"
IDL> TRIANGULATE, x, y, tri2

IDL> help

% At $MAIN$
DX FLOAT = 0.273573
DY FLOAT = 0.223977
RANGEX FLOAT = Array[2]
RANGEY FLOAT = Array[2]
TRI LONG = Array[3, 18326]
TRI2 LONG = Array[3, 19160]
X FLOAT = Array[9626]
Y FLOAT = Array[9626]
Z FLOAT = Array[9626]
 
 
Page 1 of 1    
All times are GMT
The time now is Fri Apr 18, 2014 9:03 pm