nasa = wmsfind('NASA Earth Observations','SearchField','any');
modis = refine(nasa,'land*day*month');
modis = modis(1);

server = WebMapServer(modis.ServerURL);

updatedLayer = updateLayers(server,modis);

mapRequest = WMSMapRequest(updatedLayer,server);
A = getMap(server,mapRequest.RequestURL);
R = mapRequest.RasterReference;
figure
ax = worldmap('world');
geoshow(A,R)
setm(ax,'MLabelParallel',-90,'MLabelLocation',90)
title({'MODIS Global Mosaic',modis.LayerTitle})

metadata = webread(updatedLayer.Details.MetadataURL);
disp(metadata) 

 usgsLayers = wmsfind('usgs*image','SearchField','serverurl');

serverURLs = usgsLayers.servers;
updatedLayers = [];
fprintf('Updating layer properties from %d servers.\n', ...
       numel(serverURLs));
for k=1:numel(serverURLs)
   serverLayers = refine(usgsLayers, serverURLs{k}, ...
          'SearchField','serverurl','MatchType','exact');
   serverURL = serverLayers(1).ServerURL;
   fprintf('Updating properties from server %d:\n%s\n', ...
          k, serverURL);
   server = WebMapServer(serverURL);
   try
      layers = updateLayers(server,serverLayers);
      % Grow using concatenation because layers can have any
     % length ranging from 0 to numel(serverLayers).
      updatedLayers = [updatedLayers; layers];
   catch e
       fprintf('Server %s is not responding.\n', ...
              serverURL);
       fprintf('Error message is %s\n', e.message)
    end
end