How to set map limits on mapshow?

Hi all!

I currently plotted a shapefile (attached "Expot.shp") on to the world map. However, I wanted to shorten the limits. Here is my code:

landAreas = readgeotable("landareas.shp");                   %taken from MATLAB documentation on mapshow
land = readgeotable("landareas.shp");

S = shaperead("Expot.shp") ;
N = length(S) ; 
cmap = jet(N) ;
figure
hold on
for i = 1:N
    x = S(i).X ; y = S(i).Y ; 
    x(isnan(x))=[] ; 
    y(isnan(y))=[] ;     
    patch(x,y,cmap(i,:)) ;
end
colorbar
C = [S.Corr] ;
caxis([min(C) max(C)])
hold on
geoshow(land,"FaceColor",[0.15 0.5 0.15])

I have tried using XLim and YLim but seem to get errors. Hence, I would like to know how to get the map limits to the coordinate limits:

Latitudes: 20N to 80N

Longitudes: -180 to -30W; basically the map showing part of North America and Greenland with my shapefile.

Would appreciate any help on this. Thanks!

[Walter Roberson]

axesm() allows them to be set https://www.mathworks.com/help/map/ref/axesm.html

[Keegan Carvalho]

Thank you for your reply @Walter Roberson

I did check axesm, but my shapefile is in NAD83 coordinates (Universal Transverse Mercator). I tried several projections including mercator, but i get this image:

Is there any way to code it manually? Or how would I go about this?

Answer 1

[Cris LaPierre]

I wonder why you are using shaperead instead of readgeotable for your Expot.shp file. I would create the desired map with imposed limits using axesm. I would then plot both shapefiles using geoshow.

unzip Expot.zip
land = readgeotable("landareas.shp");

axesm('pcarree','MapLonLimit',[-85 -40],"MapLatLimit",[30 80]);
geoshow(land,"FaceColor",[0.15 0.5 0.15])

data = readgeotable("Expot.shp");
N = height(data) ; 
cmap = jet(N) ;
hold on
% for loop is necessary only to set a unique color for each geopolyshape
% If you want to set the same color for all, you can remove it.
for i = 1:N
    geoshow(data(i,:),"FaceColor",cmap(i,:))
end
hold off

colorbar
clim([min(data.Corr) max(data.Corr)])

Be aware that your shapefile will have more detailed contours than what is contained in landareas.shp. Also, note that the projection used to capture Expot.shp (you say it is UTM) is different from what the default used by geoshow (Plate Carrée). Due to what a UTM map projection means, MATLAB does not/cannot create a UTM projection map of more than a single 8-by-6 depree square. See here: https://www.mathworks.com/help/map/create-utm-maps.html.

You can see if there are other projections work better: https://www.mathworks.com/help/map/summary-and-guide-to-projections.html

[Cris LaPierre]

One more update. The colors in your map do not match the colors in the colorbar. That is partly because you map is using the default colormap (parula), but your plot is using jet.

The other reason is that you are just selecting the color based on position in Expot, when you want to be selecting based on the value of Corr. This version corrects for those issues.

unzip Expot.zip
land = readgeotable("landareas.shp");

axesm('pcarree','MapLonLimit',[-85 -40],"MapLatLimit",[30 80]);
geoshow(land,"FaceColor",[0.15 0.5 0.15])

data = readgeotable("Expot.shp");
[dCorr,ia,ic] = unique(data.Corr);
N = length(dCorr);
colormap jet
cmap = jet(N);
hold on
for i = 1:length(ic)
    geoshow(data(i,:),"FaceColor",cmap(ic(i),:))
end
hold off

cb = colorbar;
cb.Label.String = "Corr";
clim([min(data.Corr) max(data.Corr)])

[Keegan Carvalho]

Thank you @Cris LaPierre

This is also perfect and yes, y aim was to plot S.Corr. But just one concern, the coordinates (latitude-longitude) don't appear on the x and y-ticks. I cannot add them in the Property Inspector as well.

Does geoshow cut them off?

[Keegan Carvalho]

Actually, there's another issue I have @Cris LaPierre.

So the Corr values which have 0s in it are actually NaNs. Is there a way to change the colours such that the polygons with 0s have no colour?

(I somehow cannot manually do it on Property Inspector)

[Cris LaPierre]

I found this statement on the documentation page for axesm:

Map axes are standard MATLAB axes with different default settings for some properties and a MATLAB structure for storing projection parameters and other data. The main differences in default settings are:

• Axes properties XGrid, YGrid, XTick, YTick are set to 'off'.
• The hold mode is 'on'.

For map axes, the lon and lat labels are called meridian and parallel labels. Turn those on and set their spacing by setting the corresponding properties.

land = readgeotable("landareas.shp");

axesm('pcarree','MapLonLimit',[-85 -40],"MapLatLimit",[30 80],...
    'MeridianLabel','on','ParallelLabel','on');
mlabel('south')
setm(gca,"PLabelLocation",10)
setm(gca,"MLabelLocation",10)

geoshow(land,"FaceColor",[0.15 0.5 0.15])

For your second question, note that setting the facecolor property to 'none' will draw the polygon but not fill it. You could therefore modify the for loop to chec the value of Corr first, and set the 'FaceColor' property accordingly.

[Keegan Carvalho]

Thanks @Cris LaPierre this works great (had wrongly set the gca positions earlier) . With regards to the using white FaceColor for Corr values of 0, I tried the following:

data = readgeotable("Expot.shp");
[dCorr,ia,ic] = unique(data.Corr);
N = length(dCorr);
colormap jet
cmap = jet(N);
hold on
for i = 1:length(ic)
    geoshow(data(i,:),"FaceColor",cmap(ic(i),:))
if data(i,:) == 0                                      %if else for Corr values = 0
geoshow(data(i,:),"FaceColor", [1 1 1])
else
geoshow(data(i,:),"FaceColor",cmap(ic(i),:))
end
end

Got a weird map followed by an error:

Operator '==' is not supported for operands of type 'table'.

Not sure if the if I got the "if-else" loop right. I guess the indexing works a different way in Matlab.

Could you please guide me with this?

[Cris LaPierre]

The if statement (not a loop) needs to check if the value of data.Corr(i) == 0. You are currently checking if all of the values in row i are 0. You also don't need the first geoshow anymore.

If there truly is no data when Corr = 0, I prefer using 'none' for the color instead of white.

unzip Expot.zip
data = readgeotable("Expot.shp");
[dCorr,ia,ic] = unique(data.Corr);
N = length(dCorr);
colormap jet
cmap = jet(N);
hold on
for i = 1:length(ic)
%     geoshow(data(i,:),"FaceColor",cmap(ic(i),:)) <----- Not needed anymore
    if data.Corr(i) == 0     % check if Corr values = 0
        geoshow(data(i,:),"FaceColor", 'none') % I prefer 'none'
    else
        geoshow(data(i,:),"FaceColor",cmap(ic(i),:))
    end
end

[Keegan Carvalho]

@Cris LaPierre this works superbly well. Thank you so much for simplifying all of this; I've understood where I went wrong now.

Very grateful for this. Thank you!!

Answer 2

[Voss]

xlim / ylim seems to work here:

% landAreas = readgeotable("landareas.shp");                   %taken from MATLAB documentation on mapshow
land = readgeotable("landareas.shp");

% S = shaperead("Expot.shp") ;
unzip('Expot.zip')
S = shaperead("./Expot/Expot.shp") ;

N = length(S) ; 
cmap = jet(N) ;
figure
hold on
for i = 1:N
    x = S(i).X ; y = S(i).Y ; 
    x(isnan(x))=[] ; 
    y(isnan(y))=[] ;     
    patch(x,y,cmap(i,:)) ;
end
colorbar
C = [S.Corr] ;
caxis([min(C) max(C)])
hold on
geoshow(land,"FaceColor",[0.15 0.5 0.15])

xlim([-180 -30])
ylim([20 80])

% this also works:
% set(gca(),'XLim',[-180 -30],'YLim',[20 80])

What errors did you get with xlim / ylim when you tried it?

[Keegan Carvalho]

Thank you @Voss this works well.

Basically I didnlt get errors, it would end up being the whole world map even after adding the Xlim and Ylim... Which was abit strange, but I think I may have missed a line of code