How do I make a contour plot?

Your data must be in a 2D grid in order for it to be contoured. There are two 2D contour graph types: 2D contour and user-defined contour.

The standard 2D contour takes a 2D grid for the values and two 1D arrays for X and Y to define the locations of the grid nodes. If you do not specify the X and Y datasets then the contour will just float between the minimum and maximum of the domain which defaults from 1 to nx and 1 to ny.

e.g.

grid = reshape(rnd(25),5,5,1);
gridx = 1//2//2.5//3//4;
gridy = 1:5;
create Viewport page_1.viewport_1
 ( XuNfirstDiagonalPoint = (10,10) %,
   XuNsecondDiagonalPoint = (40,80) %
  );
create Domain page_1.viewport_1.domain_1;
create Graph page_1.viewport_1.domain_1.graph_1
 ( XuNcolorDataGrid = "grid",
   XuNgraphType = "2DContour",
   XuNmesh = true,
   XuNxData = "gridx",
   XuNyData = "gridy"
 );
create Viewport page_1.viewport_2
 ( XuNfirstDiagonalPoint = (60,10) %,
   XuNsecondDiagonalPoint = (90,80) %
  );
create Domain page_1.viewport_2.domain_2;
create Graph page_1.viewport_2.domain_2.graph_2
 ( XuNcolorDataGrid = "grid",
   XuNgraphType = "2DContour",
   XuNmesh = true
 );

I strongly recommend specifying the location of the grid using X Data and Y Data if you can. It is much better to have the contour plot fixed in your co-ordinate space, than floating over the whole domain.

The second graph type is the user-defined contour. This graph type is drawn on an irregular grid rather than a regular grid. As the grid is irregular you must also supply a 2D grid of x values and a 2D grid of y values.

e.g.

grid = reshape(rnd(25),5,5,1);
gridx = repeaty(1:5, 5)+reshape(rnd(25),5,5,1)-.5;
gridy = repeatx(1:5, 5)+reshape(rnd(25),5,5,1)-.5;
create Viewport page_1.viewport_1
 ( XuNfirstDiagonalPoint = (10,10) %,
   XuNsecondDiagonalPoint = (40,80) %
 );
create Domain page_1.viewport_1.domain_1;
create Graph page_1.viewport_1.domain_1.graph_1
 ( XuNcolorDataGrid = "grid",
   XuNgraphType = "2DContour"
 );
create Viewport page_1.viewport_2
 ( XuNfirstDiagonalPoint = (60,10) %,
   XuNsecondDiagonalPoint = (90,80) %
 );
create Domain page_1.viewport_2.domain_2;
create Graph page_1.viewport_2.domain_2.graph_2
 ( XuNcolorDataGrid = "grid",
   XuNgraphType = "userDefinedContour",
   XuNxDataGrid = "gridx",
   XuNyDataGrid = "gridy"
 );

BUT MY DATA IS NOT IN A GRID!

If your data is found as X Y Z values then you need to convert it into a grid. If it all possible you want to avoid interpolating your data as whenever you interpolate you lose accuracy.

If you are unsure about the structure of your XYZ data, plot it first using a scatter plot:

X = list(repeaty(1:5, 5));
Y = list(repeatx(1:5, 5));
Z = rnd(25);
create Viewport page_1.viewport_1;
create Domain page_1.viewport_1.domain_1;
create Graph page_1.viewport_1.domain_1.graph_1
 ( XuNcolorData = "Z",
   XuNgraphType = "scatter",
   XuNxData = "X",
   XuNyData = "Y"
 );

You should now be able to see if your data is arranged on a regular or irregular grid. If it is then you can use various tricks to reshape your XYZ arrays into a 2D grid. I'll come back to this in a moment.

If your XYZ points are completely irregular then you will have to interpolate.

Gsharp has two standard interpolation methods bilinear and bivariate. It is not possible to say which method works best for your data. bivariate does make sure that the interpolated grid passes through the original points, but it is limited to 250 input points.

I recommend defining the location of the grid that you want to interpolate onto first. Try and choose grid nodes which resemble the locations of the original data points. For example, if you only have a few different x points you can use the function unique to pick these out:

gridx = unique(X);

or you could use the range function to create an array of values equidistant apart with the same range as your input points:

gridx = range(X,40);
gridy = range(Z,40);

Once you have defined gridx and gridy then interpolate your data:

grid = bilinear(X,Y,Z,gridx,gridy);

You can also specify regions and barriers in your interpolation - see the manual for more details. If I am using a region - I normally interpolate over the whole region and apply the region as part of the graph definition. This means that I can contour right up to the very edge of the region.

You now have your data in a grid and you can plot it as before:

X = rnd(9);
Y = rnd(9);
Z = rnd(9);
gridx = range(X,40);
gridy = range(Y,40);
grid = bivariate(X,Y,Z,gridx,gridy);
classlims = range(Z,40);
create Viewport page_1.viewport_1;
create Domain page_1.viewport_1.domain_1
 ( XuNclassData = "classlims",
   XuNclassType = "limits"
 );
create Graph page_1.viewport_1.domain_1.graph_1
 ( XuNcolorDataGrid = "grid",
   XuNgraphType = "2DContour",
   XuNxData = "gridx",
   XuNyData = "gridy"
 );
recalc;
create Domain page_1.viewport_1.domain_2
 ( XuNxMinimum = page_1.viewport_1.domain_1.XuNxMinimum,
   XuNxMaximum = page_1.viewport_1.domain_1.XuNxMaximum,
   XuNyMinimum = page_1.viewport_1.domain_1.XuNyMinimum,
   XuNyMaximum = page_1.viewport_1.domain_1.XuNyMaximum,
   XuNclassData = "classlims",
   XuNclassType = "limits"
 );
create Graph page_1.viewport_1.domain_2.graph_1
 ( XuNcolorData = "Z",
   XuNgraphType = "scatter",
   XuNxData = "X",
   XuNyData = "Y"
 );

N.B. There are two things worth noting in the above GSL.

  1. I've used recalc (he recalc command does everything that the repaint command does, apart from actually updating the canvas). This gets Gsharp to calculate the domain limits for me to use when I create the second domain. This ensures that points in my second domain correspond to the same points in my first domain.
  2. I've used classType="limits" and classData="classlims" to make sure that the both domains have the same classes and so that blue in one plot matches blue in the other.

Finally we come to the part that I put off earlier concerning reshaping a regular set of XYZ points into a 2D grid.

If you could see from the scatter plot that your XYZ forms a regular grid then we can find the locations of the x and y gridnodes using:

gridx = unique(X);
gridy = unique(Y);

If your data is arranged in the correct order:

X
 
Y
 
Z
 
1
2
3
4
5
1
2
...
 
1
1
1
1
1
2
2
...
 
?
?
?
?
?
?
?
...

Then all you need to do is reshape your Z values using reshape:

grid = reshape(Z,nx,ny,1);
# (where nx = size(gridx) and ny = size(gridy);

If your data is organised but not in this order then the simplest thing to do is grid = reshape(Z,nx,ny,1) as above. Plot the grid you have created and compare it to a scatter of the original points

(You can use the GSL code above).

You can now change your grid as follows:

To transpose: 
grid = transpose(grid);
To flip-X 
grid = slicex(grid,nx:1);
To flip-Y 
grid = slicey(grid,ny:1);

Finally, if your XYZ is missing values or you cannot work out the order then you can always use the following script. It is slower than reshape, but it should do the job.

gridx = unique(X);
gridy = unique(Y);
nx = size(gridx);
ny = size(gridy);
grid = reshape(undef,nx,ny,1);
for i=1 to size(Z)
  r = mask(1:nx,gridx=X[i]);
  c = mask(1:ny,gridy=Y[i]);
  grid[r,c] = Z[i];
endfor

If your scatter plots shows that you have an irregular grid then work out the number of rows and columns and then:

grid = reshape(Z,nx,ny,1);
gridx = reshape(X,nx,ny,1);
gridy = reshape(Y,nx,ny,1);

If you do need to transpose the grid, remember to transpose the X and Y grid as well to preserve the XYZ correspondence.

grid = slicex(grid,nx:1);
gridx = slicex(gridx,nx:1);
gridy = slicex(gridy,nx:1);