How to mask data points outside a border using geoshow?
21 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
peterhack
am 11 Nov. 2016
Beantwortet: Srishti Gaur
am 2 Aug. 2018
Hello,
I am using the presented code of this thread to mask the data points outside the border so that only data points within are displayed. Works fine so far. However, I want to use the DisplayType surface. Thus I am using geoshow to plot the data. Unfortunately it does not fill the whole map within the borders as a buffer is left.
Any ideas on this?
Here is an example:
S = shaperead('landareas', 'UseGeoCoords', true,...
'Selector',{@(name) strcmp(name,'Australia'), 'Name'});
x = linspace(min(S.Lon), max(S.Lon), 100);
y = linspace(min(S.Lat), max(S.Lat), 100);
[x,y] = meshgrid(x,y);
a = 1;
b = 9;
z = a + (b-a).*rand(100,100);
isin = inpolygon(x,y,S.Lon,S.Lat);
z2 = z;
z2(~isin) = NaN;
lnlim = [min(S.Lon) max(S.Lon)];
ltlim = [min(S.Lat) max(S.Lat)];
lt = linspace(ltlim(1), ltlim(2), 3);
ln = linspace(lnlim(1), lnlim(2), 3);
for ii = 1:2
ltbox{ii} = lt([1 2 2 1 1]'+(ii-1));
lnbox{ii} = ln([1 1 2 2 1]'+(ii-1));
end
[lnmask, ltmask] = deal(cell(2));
for ii = 1:2
for jj = 1:2
[lnmask{ii,jj}, ltmask{ii,jj}] = polybool('-', ...
lnbox{ii}, ltbox{jj}, S.Lon, S.Lat);
end
end
ltboxall = ltlim([1 2 2 1 1]);
lnboxall = lnlim([1 1 2 2 1]);
[lnmaskall, ltmaskall] = polybool('-', lnboxall, ltboxall, S.Lon, S.Lat);
figure('color','w');
worldmap('Australia');
geoshow(y, x, z2, 'DisplayType','surface')
contourcmap('jet',round(min(z2(:)),2,'significant'):1:round(max(z2(:)),...
2,'significant'),'colorbar','on','location','vertical')
for ii = 1:4
patchm(ltmask{ii}, lnmask{ii}, 'w', 'edgecolor', 'none');
end
plotm(ltmaskall, lnmaskall, 'k');
0 Kommentare
Akzeptierte Antwort
Kelly Kearney
am 15 Nov. 2016
You seem to have combined both of my suggested methods into one in this example. Are you trying to do the mask-with-NaNs option (which is easier, but may leave you with some jagged edges near the country border)? Or the patch overlay?
Either way, I think the "not filling" issue you're seeing is an artifact of how worldmap defines Australia versus how the shapefile does. The former includes a wider longitude range, probably accounting for a small island or two off the east coast somewhere (my Australian geography leaves something to be desired). If instead you call worldmap with your data limits, you should see better results:
worldmap(ltlim, lnlim);
4 Kommentare
Kelly Kearney
am 16 Nov. 2016
Bearbeitet: Kelly Kearney
am 22 Nov. 2016
Here's a patch overlay example. I decided to plot the patch in pre-projected coordinates for two reasons:
- patchm can't handle face/vertex syntax, which is the easiest way to plot patches with holes.
- This geographic region is large enough that the curvature introduced by the map projection is important. If we ignore that and try to create a mask just using straight lines in lat/lon space, then we won't completely mask the edges of the data.
% Australia border, and limits around it
S = shaperead('landareas', 'UseGeoCoords', true,...
'Selector',{@(name) strcmp(name,'Australia'), 'Name'});
ltlim = [min(S.Lat) max(S.Lat)] + [-1 1];
lnlim = [min(S.Lon) max(S.Lon)] + [-1 1];
% Some coarse fake data
nx = 20;
ny = 20;
x = linspace(lnlim(1), lnlim(2), nx);
y = linspace(ltlim(1), ltlim(2), ny);
[x,y] = meshgrid(x,y);
a = 1;
b = 9;
z = a + (b-a).*rand(ny,nx);
% Plot data as surface
figure('color','w');
worldmap('Australia');
geoshow(y, x, z, 'DisplayType','surface', ...
'zdata', ones(size(x))*-1, ... % keep below gridlines
'cdata', z, ...
'facecolor', 'flat'); % to show coarseness
geoshow(S, 'edgecolor', 'k', 'facecolor', 'none');
% Create mask overlay
ltbox = ltlim([1 2 2 1 1]);
lnbox = lnlim([1 1 2 2 1]);
[ltbox, lnbox] = interpm(ltbox, lnbox, 1);
[xbox, ybox] = mfwdtran(ltbox, lnbox); % Box around data
[xaus, yaus] = mfwdtran(S.Lat, S.Lon); % Australia polygon
[xmask, ymask] = polybool('-', xbox, ybox, xaus, yaus); % box w/ Aus. hole
[f,v] = poly2fv(xmask, ymask);
hmask = patch('faces', f, 'vertices', v, ...
'facecolor', 'w', 'edgecolor', 'none');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/185814/image.png)
Weitere Antworten (2)
peterhack
am 17 Nov. 2016
2 Kommentare
Kelly Kearney
am 18 Nov. 2016
Yeah, that's a known issue since the release of HG2 graphics in R2014b. In an attempt to be more efficient in the export of some edge case very complicated polygons, Matlab decided to go triangulation crazy; the resulting exported graphics lead to terrible aliasing artifacts in almost every pdf viewer out there.
The Mathworks doesn't consider this a bug. Everyone else who uses Matlab does. Seems to be at a stalemate.
Srishti Gaur
am 2 Aug. 2018
Dear Kelly Kearney I have a doubt, I am making a spatial plot (see image)
</matlabcentral/answers/uploaded_files/127390/Spatialplot.jpg>, Can I do something that this plot will cover my basin boundary only?
0 Kommentare
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!