Different colors for different ranges of data

figure
latlim = [-15 15];
longlim = [90 125];
worldmap(latlim, longlim)
load coast
geoshow(lat, long)

samplefactor = 1;

[Z, refvec] = etopo('etopo1_bed_c_i2.bin',samplefactor,latlim,longlim);
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z, 256);
geoshow('landareas.shp', 'FaceColor', 'none', ...
   'EdgeColor', 'black');

load Sumatra_catalog.txt

    plotm(Sumatra_catalog(:,1),Sumatra_catalog(:,2),'k.') 
    %depth = (Sumatra_catalog(:,3))

What this code does is loads a text file containing earthquake latitude, longitude and depth data. It plots the data onto a coastline as block dots.. What I want to do, is change the colour of the dots based on the different depths. For example, black for depths of 0-100, red for depths of 101-200, blue for 201-300 and so forth.

I was thinking of using a for loop with else/if statements, I've had a look at the documentation but I'm not sure how how to implement.

Any help appreciated.

Answer 1

[Walter Roberson]

Change

    plotm(Sumatra_catalog(:,1),Sumatra_catalog(:,2),'k.') 
    %depth = (Sumatra_catalog(:,3))

to

pointsize = 20;    %adjust as needed
depth = (Sumatra_catalog(:,3));

cmap = [0 0 0; %0 - 100
        1 0 0; %101-200
        0 0 1;]; %201-300
depthind = min( max(0, floor((depth-1) / 100) ), size(cmap,1)-1) + 1;   %first bin is wider than the rest
scatterm(Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, cmap(depthind,:));

[JDilla]

Hi, thanks for replying. Doesn't run, get a load of errors. Is there a simpler way using loops and if statements?

[Walter Roberson]

What error messages show up? My MATLAB session is busy at the moment so I cannot create some fake data to test with right now.

[JDilla]

Error using scatter (line 81) CData must be an RGB triplet, an M-by-1 vector of the same length as X, or an M-by-3 matrix.

Error in internal.mapgraph.GeoScatterGroup (line 70) scatter(x, y, sizedata, cdata, ...

Error in scatterm (line 62) mapgraph = internal.mapgraph.GeoScatterGroup(ax, lat, lon, args{:});

Error in practical4 (line 26) scatterm(Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, cmap(depthind,:));

>>

[Walter Roberson]

I do not see anything obvious on review, and my session is still busy. Could you post

size(depth), size(depthind), size(cmap(depthind,:))

[JDilla]

Sorry, I don't understand what you mean? I'm not a regular user of matlab, I use it sporadically so I'm not good at it. I don't understand your code, particularly

depthind = min( max(0, floor((depth-1) / 100) ), size(cmap,1)-1) + 1

[Walter Roberson]

I do not have your data, so I will create some random data:

%set up that would not apply in your actual situation
Sumatra_catalog = [rand(300,2), 500*rand(300,1)];
axesm();    %need a map axes to scatterm() into

%now the work, which _would_ apply

%extend cmap with appropriate colors
cmap = [0 0 0; %0 - 100
      1 0 0; %101-200
      0 0 1;]; %201-300
depth = Sumatra_catalog(:,3);    %will be 300 x 1
depth_bin = floor((depth - 1)/100);   %the -1 is so that 200 is in the same bin as 199
depth_bin(depth_bin < 0) = 0;    %the first bin starts at 0 instead of at 1 and because of the -1 it got converted to bin -1
depth_bin = depth_bin + 1;       %now it is a colormap index
depth_bin(depth_bin > size(cmap,1)) = size(cmap, 1);  %clamp it in case it is too large
color_triples = cmap(depth_bin, :);  %transform the index into an RGB triple

pointsize = 20;   %adjust as needed
scatterm( Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, color_triples );

[JDilla]

I've attached the data here. It follows the format (latitude, longitude, depth)

I don't understand your code, it's now bringing up a menu for 'Projection control'

Surely there's a simpler way of doing this?

[Walter Roberson]

I do not have your etopo1_bed_c_i2.bin so I cannot use your code directly.

When the bit about projection control comes up, just click on Apply. That comes up because I have not provided any parameters to axesm() . axesm() is being called on your behalf in the code that you have; the first two executable lines I show are there to create fake data and fake mapping axes to use for the real work.

[Walter Roberson]

Yes, there is an alternative code:

cmap = [0 0 0; %0 - 100
      1 0 0; %101-200
      0 0 1;]; %201-300
pointsize = 20;
scatterm( Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, Sumatra_catalog(:,3) );
caxis([0 300]);   %the upper limit should match the upper limit appropriate for the color table
colormap(cmap);

This approach will work provided that the range divides exactly equally into the number of colors in the colormap. For example if you wanted a level which was 301-350 (half the size of the others) then it would not work. The approach of creating the index, that I posted above, would be needed in such a case.

An alternate approach to creating an index would be to use histc():

cmap = [0 0 0; %0 - 100
      1 0 0; %101-200
      0 0 1;]; %201-300
[~, depth_bin] = histc( Sumatra_catalog(:,3), [0, 101, 201, 301]);
color_triples = cmap(depth_bin, :);  %transform the index into an RGB triple
pointsize = 20;   %adjust as needed
scatterm( Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, color_triples );

This has the advantage of making it easy to tweak the boundaries

[JDilla]

Instead of using a colour map, can I just define the colours myself? I don't necessarily want a range of colours, I want to determine them myself - it makes analysis easier. For example, depths of 0-100 plot as blue, depths of 101-200 plot as red. It will be easier this way, because I can then change the values of depth so maybe later I'll want 0-50 as red, and then introduce more colours as I split the data up even more. I'm looking for a pattern so I need to do it this way.

Thanks for the help.

[Walter Roberson]

cmap = [0 0 1; %0 - 100   %this is blue
      1 0 0; %101-200     %this is red
      0 0 0;]; %201-300   %this is black

define each entry in the cmap to be the exact RGB components you want for each level (double precision numbers in the range 0 to 1).

If you plan on creating new entries later that might be of different sizes then I recommend the histc() approach. Or if you are using a newer MATLAB, replace histc() with histcounts() without changing anything else. For example,

[~, depth_bin] = histc( Sumatra_catalog(:,3), [0, 50, 101, 201, 301, 489]);

cmap = [0 0 0; %0 - 50    %this is black
       1 0 0; %51-100     %this is red
       1 1 0; %101-200    %this is yellow
       0 0 1; %201-300    %this is blue
       0 1 0;]; %301-489  %this is green

[JDilla]

I still can't get this to work!

[JDilla]

is there a way of doing this with an if statement? I'm finding this way too complicated.

[Walter Roberson]

Change your code

    plotm(Sumatra_catalog(:,1),Sumatra_catalog(:,2),'k.') 
    %depth = (Sumatra_catalog(:,3))

to

[~, depth_bin] = histc( Sumatra_catalog(:,3), [0, 50, 100, 200, 300, inf]);

cmap = [0 0 0;  %0 <= depth < 50    %this is black
        1 0 0;  %50 <= depth < 100  %this is red
        1 1 0;  %100 <= depth < 200 %this is yellow
        0 0 1;  %200 <= depth < 300 %this is blue
        0 1 0;];%300 <= depth       %this is green

pointsize = 20;
scatterm(Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, cmap(depth_bin,:));

The inf in the histc is needed because you have values above 300 and they have to go into some color range. (Your largest depth is 675.1)

[JDilla]

clear
figure
latlim = [-15 15];
longlim = [90 125];
worldmap(latlim, longlim)
load coast
geoshow(lat, long)

samplefactor = 1;

[Z, refvec] = etopo('etopo1_bed_c_i2.bin',samplefactor,latlim,longlim);
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z, 256);
geoshow('landareas.shp', 'FaceColor', 'none', ...
   'EdgeColor', 'black');

load Sumatra_catalog.txt

      [~, depth_bin] = histc( Sumatra_catalog(:,3), [0, 50, 100, 200, 300, inf]);
  cmap = [0 0 0;  %0 <= depth < 50    %this is black
          1 0 0;  %50 <= depth < 100  %this is red
          1 1 0;  %100 <= depth < 200 %this is yellow
          0 0 1;  %200 <= depth < 300 %this is blue
          0 1 0;];%300 <= depth       %this is green
  pointsize = 20;
  scatterm(Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, cmap(depth_bin,:));

Error using scatter (line 81) CData must be an RGB triplet, an M-by-1 vector of the same length as X, or an M-by-3 matrix.

Error in internal.mapgraph.GeoScatterGroup (line 70) scatter(x, y, sizedata, cdata, ...

Error in scatterm (line 62) mapgraph = internal.mapgraph.GeoScatterGroup(ax, lat, lon, args{:});

Error in practical4 (line 32) scatterm(Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, cmap(depth_bin,:));

>>

[Walter Roberson]

zip together your etopo1_bed_c_i2.bin and landareas.shp and Sumatra_catalog.txt exactly as you are using, so that I can test your exact code. I have tested repeatedly on my side with the data you attached and it works here.

Also, please show the results of

size(depth_bin)
size(cmap)

[JDilla]

etopo1_bed_c_i2.bin is 300mb compressed. the landareas.shp is built in. How do I show results of size(depth_bin) ?

[Walter Roberson]

After you have run the code and had it fail, go to the command line and type in

size(depth_bin)

and press return.

300 Mb compressed is not a problem for me.

[JDilla]

 size(depth_bin)

ans =

     1     6

size(cmap)

ans =

     5     3

How should I send you the file? I cannot attach it here, the limit is 5mb.

[Walter Roberson]

You must be using an older MATLAB release. The fix is easy:

scatterm(Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, cmap(depth_bin(:),:));

[JDilla]

I get the error Subscript indices must either be real positive integers or logicals.

Error in practical4 (line 33) scatterm(Sumatra_catalog(:,1), Sumatra_catalog(:,2), pointsize, cmap(depth_bin(:),:));

I'm using matlab r2016a

[Walter Roberson]

If the above does not work,

Click on my profile and use the Contact button to send me a message so that I have your email. I will reply and you can then attach your file to the email.

Please also indicate which MATLAB version you are using so I can check with that older version.

[JDilla]

done

[Walter Roberson]

Well it turns out to be an R2016a bug in the Mapping Toolbox affecting the use of scatterm(). The bug is fixed in R2016b. To fix the bug you can

edit( fullfile( toolboxdir('map'), 'map', '+internal', '+mapgraph', 'GeoScatterGroup.m'))

and go to line 66. Where you find

                  h.LatData, h.LonData, h.SizeData, h.CData);

change that to

                  h.LatData, h.LonData, h.CData, h.SizeData);

[Walter Roberson]

If you prefer not to edit the Mathworks code, then:

clear
figure
latlim = [-15 15];
longlim = [90 125];
worldmap(latlim, longlim)
load coast
geoshow(lat, long)

samplefactor = 1;

[Z, refvec] = etopo('etopo1_bed_c_i2.bin',samplefactor,latlim,longlim);
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z, 256);
geoshow('landareas.shp', 'FaceColor', 'none', ...
   'EdgeColor', 'black');

load Sumatra_catalog.txt

[~, depth_bin] = histc( Sumatra_catalog(:,3), [0, 50, 100, 200, 300, inf]);
  cmap = [0 0 0;  %0 <= depth < 50    %this is black
          1 0 0;  %50 <= depth < 100  %this is red
          1 1 0;  %100 <= depth < 200 %this is yellow
          0 0 1;  %200 <= depth < 300 %this is blue
          0 1 0;];%300 <= depth       %this is green
pointsize = 20;

for K = 1 : size(cmap,1)
  mask = depth_bin == K;
  scatterm(Sumatra_catalog(mask,1), Sumatra_catalog(mask,2), pointsize, cmap(K,:));
end

[JDilla]

Brilliant! It worked. I'm just wondering if you could help explain what size(cmap,1) mask and depth_bin are or do?

Thank you very much for your patience and time.

[Walter Roberson]

size(cmap,1) is the number of entries in the color map.

mask is a logical vector of which entries match the current bin number. It is used for logical indexing; see http://blogs.mathworks.com/steve/2008/01/28/logical-indexing/

You can express the breakpoints of coloring your depth in an ordered sequence, and then it can make sense to talk about which case it is. Is it the "first" case, 0 <= depth < 50 ? Is it the "second" case, 50 <= depth < 100? Is it the "third" case, 100 <= depth < 200 ? And so on. So you can talk about Case #1, Case #2, Case #3, etc, with respect to each depth in your input data. The depth_bin variable records, for each depth, what the case number is. The case number of the depth is the same as the question of which row in the color table it should use: first case uses first row, second case uses second row, and so on. So the loop is going through all the different cases, picking out the entries that correspond to that case, and drawing them in the color that has been defined for that case.

[JDilla]

Thank you.

One final thing, how do I add a colourbar that corresponds to each depth?

I've tried using

colourbar
lcolorbar(labels,'fontweight','bold');

to generate it and label it but it's not working.

[Walter Roberson]

Change

[~, depth_bin] = histc( Sumatra_catalog(:,3), [0, 50, 100, 200, 300, inf]);

to

bin_edges = [0, 50, 100, 200, 300, inf];
nominal_max = 500;   %value that the data should not normally exceed
[~, depth_bin] = histc( Sumatra_catalog(:,3), bin_edges);

and after the loop that does the scatterm:

label_positions = bin_edges;
label_positions(end) = nominal_max;
caxis(bin_edges([1 end]));
h = colorbar;
h.Ticks = label_positions;

[JDilla]

Error using scribe.colorbar/set Bad value for axes property: 'YLim' Values must be increasing and non-NaN

Error in scribe.colorbar/methods>setConfiguration (line 748) set(cbar,'YLim',lim,'XLim',[0 1]);

Error in scribe.colorbar/methods (line 19) feval(args{:});

Error in scribe.colorbar (line 103) methods(h,'setConfiguration', ax);

Error in colorbar>make_colorbar (line 226) ch=scribe.colorbar(peeraxes,location,position,pvpairs{:});

Error in colorbar (line 190) [c,msg] = make_colorbar(peeraxes,location,position,pvpairs);

Error in (line 54) h = colorbar;

[Walter Roberson]

caxis(label_positions([1 end]));   %not bin_edges

[JDilla]

this has the same effect that inserting a colour bar does once the figure is generated. It makes a colour bar for the bythmetry, not the quake depths. It also changes the bythemetry files somehow and removes topography.

should I include the colour bar inside the for loop?

[Walter Roberson]

Note: at the time I wrote the above, there was no public bug report for this issue with scatterm(). The bug report has recently been made public and is https://www.mathworks.com/support/bugreports/1359012