another question has come up about hist3()

OK. So, I have access to the statistics toolbox now. And, I was wondering how I would go about plotting the hist3 function for the same code from:

http://www.mathworks.com/matlabcentral/answers/224847-how-to-create-bins-that-are-then-used-to-create-surf-function

Because [N,C]=hist3(...) doesn't actually display a plot, and having something along the lines of the Seamount:Data Point Density Histogram and Intensity plot (link is below) would be extremely helpful.

http://www.mathworks.com/help/stats/hist3.html?refresh=true

Thanks again!

Answer 1

[Star Strider]

This works, but I’m not terribly happy with the result:

[N,C] = hist3([C12 S12], {0:359; -89:89});

figure(1)
bar3(N)
grid on
xtk = linspace(-89,89, 11);
set(gca, 'XTick',xtk+90,'XTickLabel',xtk)
xlabel('DEC (-89 +89)')
ylabel('RA (1 360)')
axis tight

[jgillis16]

Sorry for the delay in response, as I was traveling. I ran the code, and I produced a plane that was uniform? It does not show the bins we imposed on the scatter plot.

[Star Strider]

I prefer to use mesh rather than bar3, although for a histogram, bar3 is more appropriate.

With this code:

fid = fopen( 'jgillis16 GWGCCatalogrm.txt');
trashLine = fgets(fid); %Skips the first line
data = textscan(fid, '%f%s%f%f%s%f%f%f%f%f%f%f%f%f%f%f%f%f', 'Delimiter', '|', 'TreatAsEmpty','~');
fclose(fid);
M1=data;
GalList.pgc = data{1};
GalList.name = data{2};
GalList.ra = data{3};
GalList.dec = data{4};
GalList.major = data{7};
GalList.abs_mag = data{14};
GalList.dist = data{15};
GalList.mass =(10.^( (-20.8 - GalList.abs_mag)./ 2.5 )) ; % Converts Absolute Magnitude into (blue luminosity) mass in units of Milky Way (luminosity) mass (how many potential forming stars in region using the blue luminosity which is rough estimate of mass).
GalList.ra = GalList.ra.*15;
GalList.dec = GalList.dec ;
GalaxyList = GalList;
C1 = GalList.ra;
S1 = GalList.dec;
C12 = round(C1);
S12 = round(S1);

[N,C] = hist3([C12 S12], {0:359; -89:89});

figure(1)
mesh(N)
grid on
xtk = linspace(-89,89, 11);
set(gca, 'XTick',xtk+90,'XTickLabel',xtk)
xlabel('DEC (-89 +89)')
ylabel('RA (1 360)')
axis tight

I get this plot:

[jgillis16]

OK. Well, I found out that the cluster I am working off of that has matlab with the stats toolbox has a graphics restriction (which means I cannot visually reproduce that plot). So, yes that was what I was looking for!

Another question: If I needed to see how many unclassified galaxies were in each bin, how would I code that? Meaning, taking each of GWGCCatalogrm.txt (classified) and NewGalaxy.txt (unclassified) right ascensions and declinations, I would be plotting them the same way, but just dividing the unclassified matrix by the classified in order to see in each bin what is the ratio of classified to unclassified. Hope that makes sense.

Thanks again!

[Star Strider]

I don’t understand ‘dividing the unclassified matrix by the classified’.

If you are doing some sort of ratio or percentage calculation, both matrices have to have the same number of elements, so you would have to use the same essential histogram calls with your two data sets. There would likely be a large number of NaN values, so you would have to decide how you want to deal with them. After that, do element-wise division (using (./)) to divide one matrix by the other to get the ratio, or the sum of the two matrices by one to get a percent.

[jgillis16]

I meant that in each square degree of the plot, I need to see how many unclassified out of classified galaxies there are. So, it would be a percentage calculation in a sense, except I would have a concrete decimal number instead. I would want all the NaN values to report as 0 (which would essentially tell me that those certain bins that display a NaN value don't have any unclassified galaxies in them).

[jgillis16]

Would this be possible? Or would I need to use a different function completely? For simplicity, I tried going back to using the surf function by combining the RA and DEC into a single matrix, then dividing that matrix by the RA and DEC of NewGalaxy as a single matrix (but since they are two different values it wouldn't work that way...obviously we would need to edit the sizes of the matrices in order for us to manipulate them) But, would this be on the right track? Or is there a simpler way?

[Star Strider]

Your data are not easy to analyse. This time, I believe it’s due to inconsistencies between your two data files and my ignorance of astronomy.

I cannot get the DEC to display correctly on these plots. I invite you to look over this code to see if you can see the problem.

The code works. The problem is likely that I don’t understand your data and how to work with it. Please be specific as to what data from what files you want to analyse.

fid = fopen( 'jgillis16 NewGalaxy.txt');
trashLine = fgets(fid); %Skips the first line
data = textscan(fid, '%f%s%f%f%s%f%f%f%f%f%f%f%f%f%f%f%f%f', 'Delimiter', '|', 'TreatAsEmpty','~');
fclose(fid);
M2=data;
C12_NG = round(M2{3});
S12_NG = round(M2{4});

fid = fopen( 'jgillis16 GWGCCatalogrm.txt');
trashLine = fgets(fid); %Skips the first line
data = textscan(fid, '%f%s%f%f%s%f%f%f%f%f%f%f%f%f%f%f%f%f', 'Delimiter', '|', 'TreatAsEmpty','~');
fclose(fid);
M1=data;
C12_GW = round(M1{3});
S12_GW = round(M1{4});

[NNG,CNG] = hist3([C12_NG S12_NG], {0:359; -89:89});    % NewGalaxy
[NGW,CGW] = hist3([C12_GW S12_GW], {0:359; -89:89});    % GWGCCatalogrm

Add = NNG + NGW;                                        % Summed Histograms
Rat = NNG ./ NGW;                                       % Ratio: NG/GW
Pct = NNG ./ Add;                                       % Percent: NG/Total

Q1 = CNG{1};                                            % Check Values
Q2 = CNG{2};                                            % Check Values

[XCtrs, YCtrs] = meshgrid(CNG{2}, CNG{1});

figure(1)
surf(XCtrs, YCtrs, Pct)
grid on
xtk = linspace(-89,89, 9);
set(gca, 'XTick',xtk,'XTickLabel',xtk)
ytk = linspace(0, 359, 10);
set(gca, 'YTick',ytk, 'YTickLabel',ytk)
qxt = get(gca, 'XTick');
qyt = get(gca, 'YTick');
xlabel('DEC (-89 +89)')
ylabel('RA (1 360)')
zlabel('Fraction:  NewGalaxy/Total')
axis tight
shading interp

Q3 = XCtrs(1:36:360, 1:18:179);                         % Check Values
Q4 = YCtrs(1:36:360, 1:18:179);                         % Check Values

NNG = accumarray([C12_NG+1, S12_NG+90], 1);
NGW = accumarray([C12_GW+1, S12_GW+90], 1);

Add = NNG + NGW;                                        % Summed Histograms
Rat = NNG ./ NGW;                                       % Ratio: NG/GW
Pct = NNG ./ Add;                                       % Percent: NG/Total

figure(2)
surf(Pct)
grid on
xtk = linspace(-89,89, 9);
set(gca, 'XTick',xtk,'XTickLabel',xtk)
ytk = linspace(0, 359, 10);
set(gca, 'YTick',ytk, 'YTickLabel',ytk)
qxt = get(gca, 'XTick');
qyt = get(gca, 'YTick');
xlabel('DEC (-89 +89)')
ylabel('RA (1 360)')
zlabel('Fraction:  NewGalaxy/Total')
axis tight
shading interp

[jgillis16]

Hey Star, thanks so much. I cannot run the data myself, as I am unable to produce the graph due to the graphical restrictions on the cluster I am working off of. I am currently waiting on someone to produce the graph for me so that I could take a look at it.

Another question. If I needed to just focus on the classified galaxies (so GWGCCatalogrm.txt) and the morphologies (so the arrays would match this time), have the same hist3 plot, and only have those galaxies' RA and DEC plotted on the x and y axis, while having the number of galaxies on the z axis, could it be done?

So my next step is to only analyze the RA and DEC (so columns 3 and 4) of GWGCCatalogrm.txt, as well as impose the morphologies of the individual galaxies (so column 5) on the hist3 plot.

[Star Strider]

My pleasure.

With respect to a histogram with the RA and DEC on the x and y axes respectively, and the counts on z, yes. That’s what I hope my latest code here does!

The morphologies will likely not be suitable for a hist3 plot, but there are others (probably stem3 or less usefully scatter3) that could do this. The problem with using hist3 is that it can only plot counts by x and y, and the others will plot whatever you want by x and y. The reason I like stem3 is that it draws stems from the plotted points to the (x,y) plane, so it’s much easier to see what goes where.

You would have to figure out how you want to designate the various morphologies, amplitude and colour being the easiest to code.

You might consider investing in the student edition of MATLAB, choosing the appropriate toolboxes. Everything we’re doing here can be done with it as far as I know, including the graphics. Might make it a bit easier on your end.

[jgillis16]

OK. Could you perhaps attach the plot you produce with the updated code provided?

I am moving the morphology question to another string: http://it.mathworks.com/matlabcentral/answers/228451-applying-a-fourth-variable-to-a-stem3-plot

[Star Strider]

I would, but since there are discrepancies between the corresponding values in the two files, it doesn’t make any sense. That’s the reason I didn’t post it here. If you look though your files to be sure the declinations and right-ascensions all span the same ranges (-89 +89 and 0 359 as appropriate), I will run my code and plot it.

I posted an Answer to your other Question. It tangentially relates to this one.