How to create bins that are then used to create surf function

I am trying to create a scatter plot that plots the right ascension against the declination and puts them into 1 degree bins. I have a file and code already that breaks down and takes the right ascension and declination into respective degree variables.

Below is the code:

load GWGCCatalogrm.txt %loads text into workspace
readCatalog( GWGCCatalogrm )
fid = fopen( '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);
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.ra = GalList.ra.*15;
GalList.dec = GalList.dec ;
GalaxyList = GalList;
C1 = GalList.ra;
S1 = GalList.dec;
c=linspace(-6,10,length(C1));

load NewGalaxy.txt
readCatalog( NewGalaxy )
fid = fopen( '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);
GalList.pgc = data{1};
GalList.name = data{2};
GalList.rad = data{3};
GalList.decd = data{4};
GalList.major = data{7};
GalList.abs_mag = data{14};
GalList.dist = data{15};
GalList.rad = GalList.rad.*15;
GalList.decd = GalList.decd ;
GalaxyList = GalList;
D1 = GalList.rad;
F1 = GalList.decd; 
s=linspace(-6,10,length(D1));
scatter(C1,S1,[],c,'filled')
hold on
scatter(D1,F1)
colorbar
xlim([0 350])
ylim([-90 90])

My files are attached. Thank you!

Answer 1

[Star Strider]

If your data are in fractional degrees, you can use round, fix or their friends to force them into integer degrees.

(I read your earlier post, but couldn’t figure out what you intend by ‘bin’. To me, that infers a histogram, but I then can’t figure out if you want the sum of the numbers of galaxies in a particular interval, and how you want to do that in the context of a 2D scatterplot.)

[jgillis16]

Ok, fixed that issue.

C12 = round(C1) 
S12 = round(S1)

From previous posts, I had to separate the galaxies with values defined in column 5 within 20 Mpc from the galaxies with no values defined in column 5 within 20 Mpc. Now, I need to make a surf plot that uses 2 matrices consisting of one matrix with the galaxies with classification (GWGCCatalogrm) and the other matrix will be without the classification with the galaxies (NewGalaxy) that is divided by one degree each 'bin' so that in each 1 degree of the surf plot I could count how many galaxies are within that 1 degree.

[jgillis16]

Side note: From the two matrices, we will need the RA and DEC so --> RA = M(:,3) DEC = M(:,4)

[Star Strider]

You are apparently categorising them by both RA and Asc, so I would need to know if you want counts that share specific RA and Asc values (so that all that have RA = 30° and Asc = 15° are summed, for instance) or counted by RA or Asc individually (RA = 30° regardless of Asc for instance) or something else?

Using the hist3 function (Statistics Toolbox) might be the easiest way:

[N,C] = hist3([C12, S12], [360, 360]);

If this doesn’t do what you want, the alternative are three nested for loops that will take at least several minutes to run.

[jgillis16]

Yes, it should be the first scenario: RA = 30° and DEC = 15° are summed.

Unfortunately, I do not have the statistics tool box. If you produce a plot, could you attach it?

And, if it's not too much trouble, could you show me how you would write the three nested for loops? I don't mind the wait time, as long as it gets the job done and I learn something so I could stop asking so many questions :)

[Star Strider]

No problem with your questions!

Use accumarray instead, and see if it gives you the output you expect:

A = accumarray([C12+1, S12+1], 1);

figure(1)
mesh(A)

It should do the same as hist3 or the nested loops.

[jgillis16]

Received an error:

Error using accumarray First input SUBS must contain positive integer subscripts.

Error in skymaskfinal (line 57) A = accumarray([C12+1, S12+1], 1);

[Star Strider]

I was hoping to avoid that by adding 1 to the matrix to counter the effect of 0°.

What are the ranges of the subscripts? It will be necessary to restrict them to a range of 1°-360° if they aren’t already. (The code worked with my test data, and produced the same plot as did hist3.)

[jgillis16]

RA: 0:360 DEC: -89:89

[jgillis16]

Attached is the plot I produced. It isn't what I expected.

[jgillis16]

This is why I wanted a surf function. So that it could account for negative values of declination too for the respective galaxies (otherwise the need to correlate right ascension and declination wouldn't make sense because RA and DEC are galactic coordinates). So, the bin size could be the Z component of the surf function, while X and Y would be the RA and DEC.

Could you write out those nested functions possibly?

[Star Strider]

O.K. You will have to add 90 to DEC for the accumarray step, then change it back in the plot. You can do this in your meshgrid call most easily, using the original values for your DEC vector.

The surfc plot doesn’t care what the actual values of the meshgrid input and output are, so long as the dimensions match with the A matrix (the accumarray output). (A meshgrid call used to be required for the surf and mesh plots but it may not be in the latest releases. You may be able to use your RA and DEC vectors directly.) You will likely have to experiment with it to get the result you want.

[jgillis16]

Ok. So I toggled with it, and it looks right. Now, I need to juxtapose the colors from the morphological types that we got from GWGCCatalogrm.txt onto the accumarray (maybe through the meshgrid call?), which were all the galaxies that had a classification for column 5. The scale for the color bar is supposed to go from -6 to 10 with -6 being red and 10 being blue. How would I impose this variable, say we called this matrix C that was a 25674x1 cell?

Could the mesh function impose C by going something along the lines of mesh(A,...,C)?

[Star Strider]

The colour definitions are easy. I would use

cm = flipud(colormap(jet(16)));

to define the colours. (You can also use colorbar with your plot.)

Depending on what you want to do, there are a few ways to go with the plot. I would use scatter3, since assigning the colours to the various points is easiest with scatter3. However, each RA-DEC point on the plot would have the same colour for all the galaxies that went into it. (If there is a 1:1 correspondence, that problem disappears.) The scatter3 plot needs vectors, not matrices, but you can do that easily enough by using meshgrid to create the matrices and then creating vectors from them as RA(:), DEC(:) and A(:). They should all match up. You would have a matching vector to define the colours, with one colour for each point, that value going from -6 to +10, but defined as cm(c) where c goes from 1 to 16, with 1 corresponding to -6 and 16 corresponding to 10.

Any element of A that has no galaxies in it (a zero count) would not plot. The best way to do this is:

A(A==0) = NaN;

[jgillis16]

'The scatter3 plot needs vectors, not matrices, but you can do that easily enough by using meshgrid to create the matrices and then creating vectors from them as RA(:), DEC(:) and A(:)'

what would be the command for this?

[jgillis16]

And, for some reason, looking at the plot now (with fresh eyes) the scaling is very off. When you said, 'O.K. You will have to add 90 to DEC for the accumarray step, then change it back in the plot. You can do this in your meshgrid call most easily, using the original values for your DEC vector.' What command would you have in mind?

[Star Strider]

I might have been clearer about the conversion of matrices to vectors. If X is a matrix, you can convert it to a linearly-addressed column vector by stating it as: X(:). So if RA and DEC are matrices created by meshgrid, you can create the vectors scatter3 requires that way. The same for the A matrix accumarray produces. They should all match. So the scatter3 call would be:

scatter3(RA(:), DEC(:), A(:), cm(c))

adding the ‘cm(c)’ column vector argument if you want to specify the colours for each dot in the scatter plot. (The order of RA and DEC might need to be reversed depending on how you specified A.)

As for the DEC vector, I was intending that you create the matrices as:

[DEC,RA] = meshgrid([-89:89], [0:359]);

So long as the dimensions agree, the plots don’t care what the actual values are.

[jgillis16]

No problem. Thanks for getting back to me. I ran the code, and now I'm coming up with the following error:

Undefined function or variable 'c'.

Error in skymaskfinal (line 107)
scatter3(C12(:), S12(:), A(:), cm(c))

Are you using a specific tool box?

[Star Strider]

No specific toolbox. The ‘cm’ matrix is the output of the colormap call I referred to earlier (corrected as):

cm = flipud(colormap(jet(17)));

and ‘c’ is a vector corresponding to the -6 to +10 values, but mapped as 1 to 17. You have to create ‘c’ by adding 7 to the column vector of integer values that are originally the -6 to +10 values.

If you want to add a colorbar, you can specify that the values go from -6 to +10 by setting 'Ticks' and 'TickLabels' arguments to it, as necessary.

[jgillis16]

OK. I keep getting errors so I am just posting what I have so far code wise. Did you manage to produce a plot?

load GWGCCatalog.txt %loads text into workspace
fid = fopen( 'GWGCCatalog.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);
A = accumarray([C12+1, S12+90], 1); 
A(A==0)=NaN;
cm = flipud(colormap(jet(17)));
mesh(A)   
[S12,C12] = meshgrid([-89:89], [-0:359]);
scatter3(C12(:), S12(:), A(:), cm(c))

[Star Strider]

What variable (the -6 to +10 values) corresponds to the colours? I don’t see it referred to specifically in this thread. I thought it was the absolute magnitude but when I looked at that the values were outside that range.

[jgillis16]

Sorry. It will be column 5 of NewGalaxy.txt, so W1.

code:

load GWGCCatalog.txt %loads text into workspace
fid = fopen( 'GWGCCatalog.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.morph = data{5};
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;

load NewGalaxy.txt %loads text into workspace
fid = fopen( '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;
GalList.pgc = data{1};
GalList.name = data{2};
GalList.ra = data{3};
GalList.dec = data{4};
GalList.morph = data{5};
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;
D1 = GalList.ra;
F1 = GalList.dec;
W1 = GalList.morph;

[Star Strider]

All I get for ‘W1’ is a cell array of tildes: '~'. This is with your posted ‘NewGalaxy.txt’ file, not the old one from your previous thread (that I still have).

What mystifies me is that 'TreatAsEmpty','~' didn’t insert NaN for that as it is supposed to do.

[jgillis16]

Ah, sorry! I am losing my mind (3 hours of sleep is not good!). We can simplify everything by keeping it simple and just concentrating on the file where we removed all the '~' column 5 galaxies --> GWGCCatalogrm.txt (I have attached). Below is the code:

load GWGCCatalogrm.txt %loads text into workspace
fid = fopen( 'GWGCCatalog.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.morph = data{5};
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;
W1 = GalList.morph;

[Star Strider]

I discovered that with some experimenting.

The problem now is that the colour vector ‘c’ is not the same size as the others in the scatter3 argument. I’m not certain how to resolve that, since there are 14435 elements of A(:) that are not NaN, and c has 25673 elements. I’m not certain the mapping would work anyway, because different values of data{5} would likely map to the same value of A. You may have to go without the colours if you want to use the histogram provided by accumarray.

My revised code is:

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);
A = accumarray([C12+1, S12+90], 1);
A = A(1:360, :);
A(A==0)=NaN;
cm = flipud(colormap(jet(17)));
c = cell2mat(cellfun(@str2num, data{5}, 'Uni',0))+7;    % Define Colour References
[S12M,C12M] = meshgrid([-89:89], [-0:359]);
figure(1)
mesh(S12M, C12M, A)
axis tight
veclens = [size(C12M(:)); size(S12M(:)); size(A(:))];   % Check Vector Lengths
figure(2)
scatter3(C12M(:), S12M(:), A(:), '.')
axis tight

Everything after the cm assignment is new or significantly revised.

[jgillis16]

It looks like the 2D plot I produced before attempting this in 3D. This is what I was looking for in terms of binning (the colors are something I have to revise from my side by extracting the color completeness from a few catalogs. It was a nice attempt and I will look further into it.)

Again, thank you so much Star! My questions are done (for now :) ).

[Star Strider]

My pleasure!

I searched through the documentation for accumarray and also hist3. Neither provides an index map from the original data to the histogram information. The nested loops would likely work if you have something else you need to do while they’re running, but that also requires that you first decide on how to handle different values of data{5} (e.g. mean, median, sum, product, something else, or simply accumulating them as part of a cell array) for each bin. I have no idea how the galaxies’ sky positions relate to data{5}. If they’re closely correlated, that might not be a problem. It’s something to think about in the interim.

I’ve reserved a specific section of my archive file for your galaxy work, with their URLs, so they all have a context!