Masking interpolation outside the boundaries of an overlay

Hello,

I am trying to find out the best way to mask an interpolation on a grid that goes outside the boundaries of an overlay.

I have 306 stations randomly located around Australia and I use the griddata function to interpolate the results between the stations. I overlay a map of Australia over the top and the interpolation causes the colouring of regions outside the borders of Australia. I would like to mask the grid squares that are outside of the boundary that I overlay. Does anyone have some ideas on how would be best to do this? I assume I would need to just manually check which grid squares fall outside the boundary and then ask MATLAB to mask these squares as white for example?

The map of Australia is right at the end of the script...it is the coast_njh It is simply an overlay. The coast_njh and the grid are both rectangular. I simply want to select the grid squares that fall outside of the map boundary and mask them with a selected colour, most likely white. I'm thinking that perhaps I have to write a script that selects chosen grid squares (such as the top left grid square being 1,1...then so on finding each square that is outside the map boundary and masking the chosen squares)

I have put my script below in case that helps? The coast_njh is the overlay that apply. It is the outline of the country of Australia

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  clear all

    load rain_3d

    load allstations_long_lat

    year = 1911:1:2010 ; year = year' ;

    month = 1:12 ; month = month' ;

%

    yr = input('Enter choice of year [1911 - 2010]> ');

%  mth = input('Enter the month number [1-12]> ');

%

    ii = find(year == yr) ;

%  ij = find(month == mth) ;

%

    x = lon_lat(:,2) ; y = lon_lat(:,3) ;  % z = rainfall anomaly  deciles 

% (1-10)

%

% calculate the long term median of Spring by first changing all values found in the chosen month that are greater than 10000 to NANs...

     vv_sep = find(rain_3d(:,49:1:148,9) > 10000) ;
     values_sep = squeeze(rain_3d(:,49:1:148,9)) ;
     values_sep(vv_sep) = nan * ones(length(vv_sep),1);
  %

% first transpose the [stations,years] to become [years,stations]
  values_sep = values_sep' ;

%

% calculate the long term median of Spring by first changing all values 

% found in the chosen month that are greater than 10000 to NANs...

     vv_oct = find(rain_3d(:,49:1:148,10) > 10000) ;
     values_oct = squeeze(rain_3d(:,49:1:148,10)) ;
     values_oct(vv_oct) = nan * ones(length(vv_oct),1);

%

% first transpose the [stations,years] to become [years,stations]
  values_oct = values_oct' ;

%

% calculate the long term median of February by first changing all values 

% found in the chosen month that are greater than 10000 to NANs...

     vv_nov = find(rain_3d(:,49:1:148,11) > 10000) ;
     values_nov = squeeze(rain_3d(:,49:1:148,11)) ;
     values_nov(vv_nov) = nan * ones(length(vv_nov),1);

%

% first transpose the [stations,years] to become [years,stations]

    values_nov = values_nov' ;

%

    values = values_sep + values_oct + values_nov ;

%

    AA = nanmedian(values) ;

%  
%

% Secondly we need to determine 'z' by finding actual rainfall totals for all stations for 1910-2010 for Spring.  Here this will be defined as 'zz'.

     tt_sep = find(rain_3d(:,49:1:148,9) > 10000) ;
     zz_sep = squeeze(rain_3d(:,49:1:148,9)) ;
     zz_sep(tt_sep) = nan * ones(length(tt_sep),1);

%

% Secondly we need to determine 'z' by finding actual rainfall totals for all stations for 1910-2010 for the Spring.  Here this will be defined as 'zz'.

     tt_oct = find(rain_3d(:,49:1:148,10) > 10000) ;
     zz_oct = squeeze(rain_3d(:,49:1:148,10)) ;
     zz_oct(tt_oct) = nan * ones(length(tt_oct),1);

%

% Secondly we need to determine 'z' by finding actual rainfall totals for all stations for 1910-2010 for Spring.  Here this will be defined as 'zz'.

     tt_nov = find(rain_3d(:,49:1:148,11) > 10000) ;
     zz_nov = squeeze(rain_3d(:,49:1:148,11)) ;
     zz_nov(tt_nov) = nan * ones(length(tt_nov),1);

%

     zz = zz_sep + zz_oct + zz_nov ;

%   

% Transpose the rainfall values for spring between 1910-2011 (inclusive)so that the matrix dimensions agree with "values" matrix.

     zz = zz' ;

%

%

     load spring_AAA

%

% calculate anomalies
%   oct_AAA = oct_AAA' ;

     AN = zz - spring_AAA ;

%

% Use prctile function to find the 10%, 20%, 30%,...,100% cut-off points for values...

    p10 = prctile(AN,10,1) ;
    p20 = prctile(AN,20,1) ;
    p30 = prctile(AN,30,1) ;
    p40 = prctile(AN,40,1) ;
    p50 = prctile(AN,50,1) ;
    p60 = prctile(AN,60,1) ;
    p70 = prctile(AN,70,1) ;
    p80 = prctile(AN,80,1) ;
    p90 = prctile(AN,90,1) ;
   p100 = prctile(AN,100,1) ;

% 

% 

    load p10_spring
    load p20_spring
    load p30_spring
    load p40_spring
    load p50_spring
    load p60_spring
    load p70_spring
    load p80_spring
    load p90_spring
    load p100_spring

%

%

all_z = NaN * ones(100,306) ;

for i = 1:30600

    if AN(i) <= p10_spring(i)

      all_z(i) = 1;

     elseif AN(i) > p10_spring(i) && AN(i) <= p20_spring(i)

      all_z(i) = 2;

     elseif AN(i) > p20_spring(i) && AN(i) <= p30_spring(i)

      all_z(i) = 3;

     elseif AN(i) > p30_spring(i) && AN(i) <= p40_spring(i)

      all_z(i) = 4;

     elseif AN(i) > p40_spring(i) && AN(i) <= p50_spring(i)

      all_z(i) = 5;

     elseif AN(i) > p50_spring(i) && AN(i) <= p60_spring(i)

      all_z(i) = 6;

     elseif AN(i) > p60_spring(i) && AN(i) <= p70_spring(i)

      all_z(i) = 7;

     elseif AN(i) > p70_spring(i) && AN(i) <= p80_spring(i)

      all_z(i) = 8;

     elseif AN(i) > p80_spring(i) && AN(i) <= p90_spring(i)

      all_z(i) = 9;

     elseif AN(i) > p90_spring(i) && AN(i) <= p100_spring(i)

      all_z(i) = 10;

    end

end

%

% transpose 'z' from [stations, time] to [time,stations]

    spring_all_z = all_z' ;

%  

    save spring_all_z 'spring_all_z'

%

    load spring_all_z

%

% Select year from spring_all_z to put in map.

    z = spring_all_z(:,ii) ;

%

% GRID RESOLUTION
% Create a grid of one degree resolution
%
% generate the regular grid domain
%

    xi = 110:0.5:155 ; yi = -45:0.5:-10 ;
    [XI,YI] = meshgrid(xi,yi) ;

%

% interpolate the data

%

    ZI = griddata(x,y,z,XI,YI);

%

% plot the results
%

    pcolor(xi,yi,ZI)
    xlabel('Longitude (^o)')
    ylabel('Latitude (^o)')

%   title([num2str(ij),num2str(ii),])

    title('Mapped Median Rainfall Deciles (mm/month) for Spring 1982')

     colormap(jet(10))
     colorbar
     caxis([0 10])

% 
%

    hold on
    coast_njh
    hold off

%

      load('MyColormaps','mycmap_median_deciles')
      set(gcf,'Colormap',mycmap_median_deciles)
      saveas(gcf, 'Spring_1982_rainfall_deciles.jpg')

[Walter Roberson]

Ger, you posted this before, and I asked several questions at the time that you never answered; and now you have deleted the previous question along with my response :(

Answer 1

[Kelly Kearney]

You'd be more likely to get responses to your questions if you post code that can actually be run (yours includes several functions and datasets that we can't access), and that only includes the relevant parts.

You can mask your data using inpolygon:

Some fake data...

 x = 110:0.5:155 ; 
 y = -45:0.5:-10 ;
 z = peaks(91); z = z(1:71,:); % some fake data
 [XI, YI] = meshgrid(x,y);

... and a coastline. Replace this with your relevant coordinates from coast_njh if you don't have the Mapping Toolbox

 Cst = shaperead('landareas.shp', ...
                'usegeocoords', true, ...
                'Selector', {@(x) strcmp(x,'Australia'),'Name'});

Now mask the data

 isin = inpolygon(XI, YI, Cst.Lon, Cst.Lat);
 z(~isin) = NaN;

And plot

 pcolor(x,y,z);
 shading flat;
 hold on;
 plot(Cst.Lon, Cst.Lat, 'k');

[Poulomi Ganguli]

What if the shape file is clipped with 3-4 shape files since I do not have any regular shape files that is directly available from MATLAB.