Contourf NaNs plotted as zeros

I couldn't find anything concerning this issue already on the forums but here is the problem I am having. The code is attempting to plot water level array as a function of x, y, and time. However, only some of the xy matrices are plotting the NaN values as white while the others are plotting the NaN values as the color of zero on the color bar. It seems odd that this is arbitrary but I'm also plotting multiple contourf arrays on the same figure so this may have something to do with it. Here is the code:

for j = 43:43%length( m.Rkm_25.wpksk );

      Fig1 = figure( 'Color', 'w', 'PaperPositionMode', 'auto', ...
      'Units', 'inches', 'Position', [ 0.01 0.01 9.49 6.24 ] );

          set( gca, 'XLim', [ 200000 300000 ], 'XTick', (200000:20000:300000), ...
              'XTickLabel', ( 200:20:300 ), ...
              'YLim', [ 80000 140000 ], 'YTick', (80000:20000:140000), ...
              'YTickLabel', ( 80:20:140 ) );

          for i = 1:length( domain );

          hval = m.Rkm_25.wpksk( j, 1 );            %index for history file
          val = m.Rkm_25.wpksk( j, 3 );            %index for map file

          Awl = s.(domain{ i }).wlvl( :, :, val );
          s.(domain{ i }).rpu = s.(domain{ i }).kfu( :, :, val );    % A domain points to removed
          s.(domain{ i }).rpv = s.(domain{ i }).kfv( :, :, val );

          Awl( s.(domain{ i }).rpu == 0 ) = NaN;                 %nan pts in A domain with no water
          Awl( s.(domain{ i }).rpv == 0 ) = NaN;

          contourf( s.(domain{ i }).xcordi, s.(domain{ i }).ycordi, Awl );
          hold on;

          end

          set( gca, 'XLim', [ 200000 300000 ], 'XTick', ( 200000:20000:300000 ), ...
              'XTickLabel', (200:20:300), ...
              'YLim', [ 80000 140000 ], 'YTick', ( 80000:20000:140000 ), ...
              'YTickLabel', ( 50:30:140 ) );
          caxis( [ -2 2 ] ); 
          colorbar
          xlabel( 'X - WA South (km)', 'fontname', 'century', 'fontsize', 16 );
          ylabel( 'Y - WA South (km)', 'fontname', 'century', 'fontsize', 16 );
          str = [ 'Water Level @ ', phase_nam{ 1, m.Rkm_25.wpksk( j, 2 )}, ...
              ' Tide ref RKM 25 low Flow', ...
              datestr( s.time.mapts( val, 1 ),'mmm-dd-yyyy HH:MM' ) ];
          title( str, 'fontname', 'Century', 'fontsize', 18 );
          set( gca, 'fontname', 'century', 'fontsize', 12 );
          axes( 'parent', Fig1, 'position', [0.57 0.70 0.25 0.20] );

%         plot( m.(mstruc_nm{ 11 }).dts2( indx.st:indx.en, 1 ), ... 
%             m.(mstruc_nm{ 11 }).dhgt2( indx.st:indx.en, 1 ), ...
%             'b-', 'linewidth', 1.5 );
% 
%         dateformat = 6;
%         datetick( 'x', dateformat );
%         set( gca, 'XLim', [ starttime endtime ], ...
%             'XTick', ( starttime:5:endtime ),'YLim',[ -2 2 ] );
% 
%         hold on ;
%         text( m.(mstruc_nm{ 11 }).dts2( hval+m.(mstruc_nm{ 11 }).indpkt-1, 1 ), ...
%             m.(mstruc_nm{ 11 }).dhgt2( hval, 1 ), 'X', 'fontname', 'century', 'fontweight','bold');
end

[curoi]

Any ideas?

[Image Analyst]

How are we supposed to know which regions are NaN or a legitimate 0? In other words, which green regions are NaN and which are zero?

[curoi]

The regions of white in the second image to the east show islands that are NaN'ed and white. The same areas in the south-east (non-inundated islands more or less) exist in the first image as well but are green for some reason even though they are represented by NaN values.

[José-Luis]

How about changing NaN's to black or some other color?

[curoi]

How would I do that? I thought NaN's were supposed to automatically be white.