How to mask data points outside a border using geoshow?

Question

peterhack 2016-11-11

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/311766-how-to-mask-data-points-outside-a-border-using-geoshow

回答： Srishti Gaur 2018-8-2

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Kelly Kearney 2016-11-15

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/311766-how-to-mask-data-points-outside-a-border-using-geoshow#answer_243305

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

peterhack 2016-11-15

编辑：peterhack 2016-11-15

在 MATLAB Online 中打开

Yes, you are right that I have been mixing up the methods. However, I doubt this is a worldmap vs. shapefile problem. I have been trying both for various countries.

I will now concentrate on the patching method to get a clean result. When I use your example, it does not patch overlay all the area.

Do you have any suggestions on that?

file = 'TM_WORLD_BORDERS-0.3.shp';
S = shaperead(file, 'usegeocoords', true, 'selector', ...
    {@(x) strcmp(x,'India'), 'NAME'});
x = linspace(min(S.Lon), max(S.Lon), 100);
y = linspace(min(S.Lat), max(S.Lat), 100);
[x,y] = meshgrid(x,y);
z = rand(100);
lnlim = [65 100];
ltlim = [5 40];
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(ltlim, lnlim);
geoshow(y, x, z, 'DisplayType','surface')
for ii = 1:4
    patchm(ltmask{ii},  lnmask{ii}, 'w', 'edgecolor', 'none');
end
plotm(ltmaskall, lnmaskall, 'k');

Kelly Kearney 2016-11-15

在 MATLAB Online 中打开

Ah, found the problem. You're using a surface to plot the data:

geoshow(y, x, z, 'DisplayType','surface')

Even though it looks 2D in this view from above, the data is actually plotted as a 3D surface, with z-values >0 in this case, while all 2D objects (such as the masking patch) is set at z = 0.

If you use

geoshow(y, x, z, 'DisplayType','texturemap')

it will create a 2D raster at z=0, and the example will work.

Keep in mind that one difference between the 'surface' and 'texturemap' options is their treatment of NaNs; 'surface' plots those NaN-regions as transparent (like surface, pcolor, etc.) while 'texturemap' sets them to the minimum colorbar value (like image, imagesc, etc). So for the inpolygon NaN-masking method, use the former, but for the patch overlay method, use the latter.

I'll restate that the patch overlay is really only necessary if your trying to show coarser data, where jagged masking along the border of your polygon is really problematic. Assuming the data in your example matches your real data, I don't think that's much of a problem:

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;
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')
geoshow(S, 'edgecolor' ,'k', 'facecolor', 'none')

peterhack 2016-11-16

Thanks for all the background information, now I got the difference.

However, I already got the latest output map that you have posted with the code of my initial past. And the reason for this thread is that with the surface typ I get an ugly result with a lot of boardbuffer where just nothin is displayed. That it was bothers me and that is why I have been asking whether it is possible to combine the surface type with the patch overlay method to avoid jagged masking along the border.

The reason I want to use the surface type is that I have quite coarse data not looking great using the texturemap style. Hopefully my intention got clearer now.

So is it possible combining surface and the patch overlay method? I had a look at fill3m but do not know how to utilize this...

Kelly Kearney 2016-11-16

编辑：Kelly Kearney 2016-11-22

在 MATLAB Online 中打开

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');

请先登录，再进行评论。

Answer 2

peterhack 2016-11-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/311766-how-to-mask-data-points-outside-a-border-using-geoshow#answer_243660

Great work, thank you so much! Only issue left is that printing it as a native chart (.eps, .pdf) leaves a lot of stripes around the borders. I do not now whether this is due to the polygon or to some printing properties?

2 个评论
显示无隐藏无

Kelly Kearney 2016-11-18

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.

Read lots of related discussion, and some suggested workarounds, here.

peterhack 2016-11-18

Ah ok, seems to work out using opengl :-)

请先登录，再进行评论。

Answer 3

Srishti Gaur 2018-8-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/311766-how-to-mask-data-points-outside-a-border-using-geoshow#answer_331280

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

How to mask data points outside a border using geoshow?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（2 个）

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

How to mask data points outside a border using geoshow?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（2 个）

2 个评论 显示 无隐藏 无

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论