How can I plot overlay black land?

At the moment I have a plot that looks like the attached, with zero values being land.

How can I make the land black in this plot? without affecting the zero values in the data..

thanks, Michael

[Michael]

the data is on a 1682 x 3600 grid...

Answer 1

[Kelly Kearney]

Plot using pcolor (not image or similar), mask out the land with NaNs, and set the axis color to black.

If you don't have a specific land mask already defined, you can use inpolygon with a coastline dataset. For example

lat = linspace(-90, 90, 100);
lon = linspace(-180, 180, 100);
z = peaks(100);

[long, latg] = meshgrid(lon, lat);

C = load('coast');
island = inpolygon(long, latg, C.long, C.lat);

z(island) = NaN;
pcolor(long, latg, z);
set(gca, 'color', 'k');
shading flat;

The inpolygon call will be a bit slow with a matrix your size, but you should only have to call that once and then can save the mask for later use. I use my inpolygons function instead, and it tends to be faster for tasks like this.

[Michael]

I tried pcolor but it plots a completely black image with no data shown...the land values are now NaN's though

[Michael]

Is there a way to plot NaN's as black?

[Kelly Kearney]

Did you add the

 shading flat

command? Without it, at your resolution, the black outlines of each grid cell (which are present under the default faceted shading) will overwhelm the color. If you run my example completely (and adapt it to your data), you should see the land-NaNs plotted in black.

[Michael]

I tried changing the size of lat & lon to be the same size as my grid but when I tried to you inpolygon it took too long. Should I do this? or should I use my own mask?

[Image Analyst]

Michael said "values in the data are also black" so no amount of plotting options and trickery will solve that. As I said in my comment, Michael, you have to change your data so that there are no "land" pixels in what is supposed to be the "water region". As long as they have the same intensity value, they will be there (unless you shrink the image so much they don't show up) and no amount of plotting gymnastics and chicanery will change that.

[Kelly Kearney]

Using the method I showed above, where land values are designated as such based on coastline polygons, it shouldn't matter. I assume he means that the land values are given a default value if 0, which also happens to be a perfectly valid data values. Therefore assigning the color black to 0 is not a good solution. Calculating the real land mask (based on polygons, not data == 0) will solve the problem.

Michael, did you try the inpolygons function I linked to? It's still slow, but not as slow as inoolygon. Of course, if you already have a land mask associated with the data, of course use that.

Excuse the typos, typing on a phone and it doesn't play well with this forum at all.

[Kelly Kearney]

Alternatively, assuming you're not planning to zoom in to full resolution with this figure, there's really no reason that your land mask has to be calculated at such a high resolution.

The following example calculates the land mask at 1/4 the resolution, then upsamples it for use as a mask. Only takes a ~30 sec on my computer.

Note that I had to add in the flatearthpoly call to close off Antarctica properly. I'm assuming you have the Mapping Toolbox, since you didn't comment on my use of the coast.mat dataset.

nlat = 1682;
nlon = 3600;

lat = linspace(-90, 90, nlat);
lon = linspace(-180, 180, nlon);
z = peaks(max(nlat,nlon));
z = z(1:nlat, 1:nlon);

[long, latg] = meshgrid(lon, lat);

C = load('coast');
[latc, lonc] = flatearthpoly(C.lat, C.long);

lat2 = linspace(-90,90,nlat/4);
lon2 = linspace(-180,180,nlon/4);
[long2, latg2] = meshgrid(lon2, lat2);

island = inpolygons(long2, latg2, lonc, latc);
island = interp2(long2, latg2, island, long, latg, 'nearest');

z(island) = NaN;
pcolor(long, latg, z);
set(gca, 'color', 'k');
shading flat;
hold on;
plot(lonc, latc, 'w');
colorbar;

Answer 2

[Joseph Cheng]

You can change the colormap used to display the image. First obtain the colormap used by:

 map = colormap;

then change the color map corresponding to zero to black [0 0 0]

 map(1,:) = 0;

then apply the colormap:

 colormap(map);

[Joseph Cheng]

example:

images = 0:.01:.5;
images=repmat(images,51,1);
figure(1),imagesc(images)
map = colormap;
map(1,:)=0
figure(2),imagesc(images),colormap(map)

[Michael]

Thanks, although I have tried this previously, but when I used it for another plot where there are zero values in the ocean data as well it looks bad. Also, I don't like that it shows black = 0 in the colorbar.

[Joseph Cheng]

If that is a case you may want to check the imoverlay() function at the file exchange. using that function you would need to create a binary mask of the land values and then apply it using the imoverlay() function. I'm not sure what this will do to the color bar though.

[Michael]

I tried this function but I cannot seem to get it to work...

I am using imagesc and this example seems to be for the imshow command.

[Joseph Cheng]

well you could do something like this.

%create dummy image
images = 0:.01:1;
images=repmat(images,101,1);

%create mask in the center.
mask = images==.5;
%create a foreground flat image offset it to a value not used by the image
maskImage = max(images(:))*ones(size(images))+.1;

f1 = figure;
h1= imagesc(images), colorbar;

%generate the colormap
map = colormap;
%append black to the end of the colormap
map = colormap([map;gray(1)]);

hold on;
%overlay the earth image with the foreground mask
 h2 = imagesc(maskImage);
%using the logic mask set the foreground to be transparent.
set(h2,'alphaData',mask)

%edit the colorbar to only show values for the earth image.
cbar_handle = findobj(f1,'tag','Colorbar')
pause(.1)
set(cbar_handle,'Ylim',[0 1])

the only problem i have with this right now is that the last line needs to have a pause of some value for it to work in a script. so here i have created a mask for the value of .5. in your case i'd create a mask of the land values of zero where i know i'll only have zero values for land then use that for the other data sets where there are zero values in the sea.

Answer 3

[Image Analyst]

If you don't like the black line in the colorbar, you can convert to rgb, like this untested code:

% Create standard colormap
customColorMap = jet(256); 
% Make bottom color be black
customColorMap(1,:) = [0,0,0];
% Convert to RGB
rgbImage = ind2rgb(indexedImage, customColorMap);
% Display
imshow(rgbImage);
% Show color map, standard one, not custom one
colormap(jet(256));
colorbar;

[Michael]

thanks, although I have to make it so that just the land is black, at the moment values in the data are also black.

[Image Analyst]

You must have some pixels in the water that are the same intensity as the land. How do you know they are not land? Maybe they are islands? If you want to filter out small specks that you don't want to be black, then use bwareaopen.

I'm going to be traveling for the next week after tomorrow, and not online very much, so if you can upload your indexed image (not the color one that you uploaded already), and point out what areas you don't want black, then I'll see if I can get time to look at it.

[Image Analyst]

Did you see the comment about uploading your image (not a screenshot)?

Answer 4

[corvianawatie]

if your data is Z(i,j), where i=1682 and j=3600 then you must have land masking data (contain of 0 and 1), e.g. if land then the value is one and if ocean then the value is zero. Then do this calculation:

for i=1:1682;
for j=1:3600;
    if (land(i,j)==1);
    Z(i,j) = NaN;
    else;
    Z(i,j) = Z(i,j);
    end;
end;
end;
end;

pcolor (Z);

[Image Analyst]

I'm not sure Michael has two arrays, a land mask and an indexed image. If he did, he could do what you say in one line:

Z(land)=nan;

and I don't know why you use pcolor instead of imshow(), imagesc(), or image()? Any reason? Did you know that pcolor() does not show the last row and column, though with thousands of rows and columns, and shrinking the image for display, you'll lose some to subsampling for display so you won't notice it. But I never use pcolor() because of that reason (it chops off the last row and column).

Anyway, I think it's evident Michael is just ignoring this thread now.

[Michael]

I am not ignoring this thread at all! thanks for the suggestions. I have/am trying to use the advice given here but so far it is not working.

I have a landmask already, but how can I plot overlay keeping the original data. When I use imshow it produces the land map but deletes the data..

I currently use 'hold on', is this incorrect?

[Michael]

the land is already = NaN values now in the original data, so I guess this is the same as using a landmask?

[Michael]

I finally managed to do it (using an answer from here: http://www.mathworks.de/matlabcentral/answers/73492-assigning-different-color-to-nan-values-in-2d-matrix)

 imagesc(data) ;
 [r,c] = find(isnan(data)) ;
 hold on ; 
 plot(c, r, 'ro') ;

Thanks for all your input guys!

[Michael]

[Image Analyst]

That's a screenshot again. We don't want axes labels, tick marks, colorbars, etc. - just the image itself. Save your image variable in a .mat file and attach the .mat file.

And I don't see how "plot(c, r, 'ro')", which plots a red circle over every nan value, turns the land from dark blue to black.

[Michael]

I tried but the .mat/.fig file is too large to attach to this forum.

I interactively selected the red circles within the plot and changed the line width from 6 -> 0.5 and changed the color to black and it gave me the picture attached above.

[Kelly Kearney]

If the NaN mask is in place, I'll still argue that it's much more straightforward to just use pcolor than trying to replicate the same appearance with contortions of image and plot.

figure;
axes('color', 'k');
pcolor(data);
shading flat;

Answer 5

[Chad Greene]

I turned Kelly's inpolygon solution into a function called landmask. It's worth noting that downsampling the coast line does make computation faster, but it introduces an error: Each land mass is separated by a row of NaNs in the coast line arrays. If you skip over a NaN row, land masses get conjoined. I included a workaround to that problem in the landmask function.