Discussion:
Help creating a circular contour plot
(too old to reply)
Kate
2009-03-17 13:21:01 UTC
Permalink
Hello everyone,

I would really appreciate some help with this. This problem came up while writing a program to calculate the luminance distribution over the sky at a given time for my 4th year project (I’m a mechanical engineering student at Manchester, UK).

I have calculated the relative luminance for a sky, as modelled as a hemisphere. So I defined two angles: x=linspace(0,2*pi,360) and y=linspace(0,pi/2,360) with x being the angle of azimuth of a point and y being the angle of elevation from the horizon up to the zenith. I calculated my relative luminance for each combination in a for loop, so I have RelLum as a 360x360 matrix. What I want to do is draw a circular contour map of the Relative luminance, with the angle of elevation gong from 0 to 2pi around the edge of the plot. If you google "Typical clear sky luminance pattern in the plan projection of the sky vault normalised to zenith", the first image you see will give you an example of what I want. At the moment, my contour plot is square, with the x being on one axis and y on the other.

As it stands the data is in a square, Cartesian form and I want it to become circular as shown in the example. As RelLum is defined by a full 2pi of angles, this should be easy but I cannot figure out how to do it. If someone could help me I would be extremely grateful. I’m not an experienced Matlab users and so not an experienced poster, please let me know if there’s anything else you need from me to figure it out.
us
2009-03-17 13:35:04 UTC
Permalink
"Kate"
As it stands the data is in a square, Cartesian form and I want it to become circular as shown in the example. As RelLum is defined by a full 2pi of angles, this should be easy but I cannot figure out how to do it...
did you look at these

help cart2pol;
help polar;

us
Peter Boettcher
2009-03-17 13:58:45 UTC
Permalink
Post by Kate
Hello everyone,
I would really appreciate some help with this. This problem came up
while writing a program to calculate the luminance distribution over
the sky at a given time for my 4th year project (I’m a
mechanical engineering student at Manchester, UK).
I have calculated the relative luminance for a sky, as modelled as a
hemisphere. So I defined two angles: x=linspace(0,2*pi,360) and
y=linspace(0,pi/2,360) with x being the angle of azimuth of a point
and y being the angle of elevation from the horizon up to the
zenith. I calculated my relative luminance for each combination in a
for loop, so I have RelLum as a 360x360 matrix. What I want to do is
draw a circular contour map of the Relative luminance, with the angle
of elevation gong from 0 to 2pi around the edge of the plot. If you
google "Typical clear sky luminance pattern in the plan projection of
the sky vault normalised to zenith", the first image you see will give
you an example of what I want. At the moment, my contour plot is
square, with the x being on one axis and y on the other.
Why not provide the full link to the image?
Post by Kate
As it stands the data is in a square, Cartesian form and I want it to
become circular as shown in the example. As RelLum is defined by a
full 2pi of angles, this should be easy but I cannot figure out how to
do it. If someone could help me I would be extremely
grateful. I’m not an experienced Matlab users and so not an
experienced poster, please let me know if there’s anything else
you need from me to figure it out.
This is not particularly simple.

I would start by defining the x/y grid in cartesian space (it's the
space you'll be plotting on, so it will provide even resolution for the
plot). Then, convert the x/y pairs on that grid to az/el polar
coordinates. Compute your luminance on those coordinates. When you
contour-plot that, you'll get an answer CLOSE to your example image.

% x y defined on a -1 to 1 grid

[x y] = meshgrid(linspace(-1,1,360));
az = atan2(y, x); % radians
el = pi/2*(1-sqrt(x.^2 + y.^2)); % radians, pi/2 at center, 0 at edge

RelLum = some computation on az and el

contour(RelLum)


Problems: No labels, and an output that ranges below the horizon.

First eliminate the data that drops below the horizon (el < 0):

RelLum(el < 0) = NaN;
contour(RelLum)


OK, that's better.

I'm afraid the labels and axes you need to draw yourself. 'type polar'
will give an example of a MATLAB function that does a full
do-it-yourself polar plot. You'll either have to dig into it yourself,
or ask someone with more time than I have to pull out the necessary bits
and modify it to your display style.

Maybe there's something on the file exchange?

-Peter
Kate
2009-03-17 15:08:01 UTC
Permalink
Hi Peter,

Thanks very much for your help, I shall give it a go this afternoon.

Well pointed out that I should have put in the link, just quickly for future reference the example image can be found at Loading Image...

Kate
John D'Errico
2009-03-17 15:09:01 UTC
Permalink
Post by Kate
Hello everyone,
I would really appreciate some help with this. This problem came up while writing a program to calculate the luminance distribution over the sky at a given time for my 4th year project (I&#8217;m a mechanical engineering student at Manchester, UK).
I have calculated the relative luminance for a sky, as modelled as a hemisphere. So I defined two angles: x=linspace(0,2*pi,360) and y=linspace(0,pi/2,360) with x being the angle of azimuth of a point and y being the angle of elevation from the horizon up to the zenith. I calculated my relative luminance for each combination in a for loop, so I have RelLum as a 360x360 matrix. What I want to do is draw a circular contour map of the Relative luminance, with the angle of elevation gong from 0 to 2pi around the edge of the plot. If you google "Typical clear sky luminance pattern in the plan projection of the sky vault normalised to zenith", the first image you see will give you an example of what I want. At the moment, my contour plot is square, with the x being on one axis and y on the other.
As it stands the data is in a square, Cartesian form and I want it to become circular as shown in the example. As RelLum is defined by a full 2pi of angles, this should be easy but I cannot figure out how to do it. If someone could help me I would be extremely grateful. I&#8217;m not an experienced Matlab users and so not an experienced poster, please let me know if there&#8217;s anything else you need from me to figure it out.
You have two angles, that define a point on the
surface of a unit sphere, actually, the surface of
a unit hemisphere, since one of those angles goes
only from 0 to pi/2.

On the surface of that sphere, you have some other
function defined, for which you wish to draw the
contours. You think of it as a polar coordinates
problem, since your plot is actually a projection
of the hemisphere into the plane, and is thus a
circle.

The projection you indicate has you standing at
the center of the sphere and looking around. Here
x maps to a longitudinal polar angle, and y is the
angle above the horizon. S when y = pi/2, this
maps to the point directly overhead. To view this
in polar coordinates, you need to think of the polar
radius as

r = pi/2 - y

See that when y = pi/2, i.e., the point directly
above you, this maps to the origin (r = 0) in polar
coordinates. At the other end of the scale, when
y = 0, this results in a point on the circumference
of the circle with radius pi/2.

How do we encode this into Matlab? First, I'll
suggest that you have built too fine of a grid. 30
steps or so is probably a bit more manageable,
and will produce entirely reasonable looking plots.

n = 30;
x=linspace(0,2*pi,n);
y=linspace(0,pi/2,n);
[theta,r] = ndgrid(x,pi/2 - y);

% since matlab plots in cartesian coordinates,
% we need to convert this polar form into the
% cartesian coordinates
[xc,yc] = pol2cart(theta,r);

% Now we need to triangulate the polar domain.
[I1,I2] = ndgrid(1:n-1,1:n-1);
ind = I1(:) + (I2(:)-1)*n;
tri = [[ind,ind+1,ind+n+1];[ind,ind+n,ind+n+1]];

% We could now use a tool like trimesh or trisurf.
% For example, here is the polar mesh:
trimesh(tri,xc,yc)
axis square

% Now, plot the surface itself using trimesh, or trisurf.
trisurf(tri,xc,yc,RelLum)
axis square

For the contours, you will need to download
tricontour from the file exchange.

http://www.mathworks.com/matlabcentral/fileexchange/10408

I've not done the extra work to resolve the wrapping
problem across the 0-360 degree boundary, nor have
I coalesced the points at the polar origin into a
single point. This is all easily enough done, but not
necessary in your specific problem.

John
Kate
2009-03-17 23:43:01 UTC
Permalink
Hi John,

Thanks so much for your help. You and Peter giving me your time has made me feel all happy about the world and the good of humankind - bit soppy I know!

The trisurf command gives a very cool image. I know very little about triangulation, just that it is used in reverse engineering processes to mesh a cloud of points out, why is it necessary in this case? What does it bring to the table, as it were? I'm intrigued. Because I don't really understand it, I wasn't able to implicate the tricontour file successfully as I don't understand what it is asking for in the inputs. I've copied below what you get by typing "help tricontour" in matlab, which explains what the input variables are, could you tell me how to adapt it for this specific case? I mean, what variables do I put in for p,t,F and N? Thank you very very much.

help tricontour
Contouring for functions defined on triangular meshes

TRICONTOUR(p,t,F,N)

Draws contours of the surface F, where F is defined on the triangulation
[p,t]. These inputs define the xy co-ordinates of the points and their
connectivity:

P = [x1,y1; x2,y2; etc], - xy co-ordinates of nodes in the
triangulation
T = [n11,n12,n13; n21,n23,n23; etc] - node numbers in each triangle

The last input N defines the contouring levels. There are several
options:

N scalar - N number of equally spaced contours will be drawn
N vector - Draws contours at the levels specified in N

A special call with a two element N where both elements are equal draws a
single contour at that level.

[C,H] = TRICONTOUR(...)

This syntax can be used to pass the contour matrix C and the contour
handels H to clabel by adding clabel(c,h) or clabel(c) after the call to
TRICONTOUR.

TRICONTOUR can also return 3D contours similar to CONTOUR3 by adding
view(3) after the call to TRICONTOUR.

Type "contourdemo" for some examples.

See also, contour, clabel

Loading...