Map with Latitudes and Longitudes, filled with shades

Fun maps...

 

%%

clear all

 

% Sample latitude and longitude data

ilat = 25*rand(1,50)+25;

ilon = 50*rand(1,50)-120;

% Sample data for color-coding the points

data = 10*rand(1,50);

 

close all

latlim=[min(ilat) ,max(ilat)];

lonlim=[min(ilon) , max(ilon)];

% find number of stations lat lon and number of data

h1 = figure(1);

ax = worldmap(latlim, lonlim);

states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);

geoshow(ax, states, 'Facecolor', [1 1 0.9]);

hold on

map = colormap;

colorvector = map(1:64,:);

pmvect = linspace(0, max(data), 64);

for jj = 1:length(data)

    %jj

colorb = interp1(pmvect, colorvector, data(jj));

plotm(ilat(jj),ilon(jj), data(jj),'o', 'MarkerEdgeColor','k','MarkerFaceColor', colorb, 'MarkerSize', 5);

hold on   

end

caxis([0  max(data)])

colorbar

hold off

set(gcf,'Color','white')

Plotting NPP hdf data

A sample code to plot the NPP data

The data source:

http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php 

 filewant = 'cafe.2007001.hdf';

 

info = hdfinfo(filewant);

 

% This is how you can investigate the file structure...

info.SDS

%% note the names

info.SDS.Attributes.Name

%% note the values of the attributes

info.SDS.Attributes.Value

%% read the npp variable

npp  = double(hdfread(filewant,  'npp' ));

 

npp(npp==-9999) = NaN;

% figure;

% imagesc(npp)

npp = flipud(npp); % because npp is flipped upside down wrt lat/lon for map

 

%%

iLat = linspace(-90,90,1080);

iLon = linspace(-180,180,2160);

[Lat, Lon] = ndgrid(iLat,iLon);

 

figure

load coastlines % coastineData.mat

latlim=[floor(min(min(Lat))),ceil(max(max(Lat))) ]; 

lonlim=[floor(min(min(Lon))),ceil(max(max(Lon))) ];

ax = worldmap(latlim, lonlim);

surfacem(Lat, Lon, npp);

geoshow(coastlat, coastlon,'Color', 'k' )

colormap; set(gcf,'Color','white')

colormap jet

map2 = colormap; 

map2( 1, : ) = 1; colormap(map2);

colorbar 

 

title('global dataset')

  

 

disp('Next we will try to crop out the area of interest')

%% This is an example of how to crop out to the area of interest...

iLat2 = linspace(-15,45,1080);

iLon2 = linspace(-90,0,2160);

[Lat2, Lon2] = ndgrid(iLat2,iLon2);

 

F = griddedInterpolant(Lat,Lon,npp); % read doc griddedInterpolant

npp_crop = F(Lat2,Lon2);

 

%

 

figure

load coastlines % coastineData.mat

latlim=[floor(min(min(Lat2))),ceil(max(max(Lat2))) ]; 

lonlim=[floor(min(min(Lon2))),ceil(max(max(Lon2))) ];

ax = worldmap(latlim, lonlim);

surfacem(Lat2, Lon2, npp_crop);

geoshow(coastlat, coastlon,'Color', 'k' )

colormap; set(gcf,'Color','white')

colormap jet

map2 = colormap; 

map2( 1, : ) = 1; colormap(map2);

colorbar 

title('NPP (mgC m-2 day-1)')

Alternative code that I have seen:

 https://hdfeos.org/zoo/OTHER/npp_2010361_hdf.m

 

Create Clickable imagesc for 2D data

Looking at the figure is fun, specially if you can click and get the values. So, created this little function which can give you title and data tips @cursor location.

function imagesca(x, titlewant)
h = figure;
imagesc(x);
colorbar
try
    title(titlewant);
catch ME
    title('Figure')
end
datacursormode on

Map with inset

So, I had to zoom into NYC area with the counts of data shown for lat/lon points
here is a figure:

NYC Zoomed in to show more details

close all
latlim = [40, 45];
lonlim = [-80,-72];
% latlim = [40, 42];
% lonlim = [-80,-73];
% find number of stations lat lon and number of data
h1 = figure(1);

ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'Facecolor', [1 1 0.9]);
hold on
map = colormap;
colorvector = map(1:64,:);
pmvect = linspace(0, max(saveI), 64);
for jj = 1:41
    jj
   
    %      textm(st_lat(jj),st_lon(jj), num2str(jj));
    colorb = interp1(pmvect, colorvector, saveI(jj));
%     plotm(Lat1(jj),Lon1(jj), jj,'o', 'MarkerEdgeColor','k','MarkerSize',saveI(jj)+1);
plotm(Lat1(jj),Lon1(jj), saveI(jj),'o', 'MarkerEdgeColor','k','MarkerFaceColor', colorb, 'MarkerSize', 5);
hold on
%     textm(Lat1(jj),Lon1(jj),sprintf('%.1d', jj))
   
end
caxis([0  max(saveI)])
% colorbar
title('The number of Episodic (PM2.5>35\mug/m^{-3}) events in NY')
hold off


h2 = figure(2);
latlim = [40.5, 40.9];
lonlim = [-74.3,-73.5];
ax = worldmap(latlim, lonlim);
set(ax, 'FontSize', 1)

states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'Facecolor', [1 1 0.9]);
hold on
map = colormap;

colorvector = map(1:64,:);
pmvect = linspace(0, max(saveI), 64);
for jj = 1:41
    jj
    %      textm(st_lat(jj),st_lon(jj), num2str(jj));
    colorb = interp1(pmvect, colorvector, saveI(jj));
%     plotm(Lat1(jj),Lon1(jj), jj,'o', 'MarkerEdgeColor','k','MarkerSize',saveI(jj)+1);
plotm(Lat1(jj),Lon1(jj), saveI(jj),'o', 'MarkerEdgeColor','k','MarkerFaceColor', colorb, 'MarkerSize', 5);
hold on
%     textm(Lat1(jj),Lon1(jj),sprintf('%.1d', jj))
   
end
caxis([0  max(saveI)])
% set(findobj(gcf, 'type','axes'), 'Visible','off')
hold off
% trying to inset the figure
inset_size=0.4;
inset_size=inset_size*.7;
figure;
new_fig=gcf;
main_fig = findobj(h1,'Type','axes');
h_main = copyobj(main_fig,new_fig);
set(h_main,'Position',get(main_fig,'Position'))
colorbar
inset_fig = findobj(h2,'Type','axes');
h_inset = copyobj(inset_fig,new_fig);
ax=get(main_fig,'Position');
% set(h_inset,'Position', [.7*ax(1)+ax(3)-inset_size .5*ax(2)+ax(4)-inset_size inset_size inset_size])
set(h_inset,'Position', [0.3*ax(1) 1.5*ax(2) inset_size*2 inset_size])