How to set map limits on mapshow?

26 次查看(过去 30 天)
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!
  2 个评论
Keegan Carvalho
Keegan Carvalho 2022-6-27
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?

请先登录,再进行评论。

采纳的回答

Cris LaPierre
Cris LaPierre 2022-6-27
编辑:Cris LaPierre 2022-6-27
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
  7 个评论
Cris LaPierre
Cris LaPierre 2022-6-28
编辑:Cris LaPierre 2022-6-28
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
Keegan Carvalho 2022-6-28
@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!!

请先登录,再进行评论。

更多回答(1 个)

Voss
Voss 2022-6-27
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?
  1 个评论
Keegan Carvalho
Keegan Carvalho 2022-6-28
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

请先登录,再进行评论。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by