Documentation

Search the WMS Database for all layers with 'nasa' in the ServerURL field.

nasa = wmsfind('nasa','SearchField','serverurl');

Use the WMSLayer.refine method to refine your search to include only those layers with the phrase 'bluemarbleng' in the LayerName field. This syntax creates an exact search.

layer = nasa.refine('bluemarbleng', 'SearchField','layername', ...
   'MatchType','exact');

Use the wmsread function to retrieve the Blue Marble Next Generation layer.

[A,R] = wmsread(layer);

The wmsread function returns A, a geographically referenced raster map, and R, a raster reference object that ties A to the EPSG:4326 geographic coordinate system. The geographic limits of A span the full latitude and longitude extent of layer.

Open a figure window, set up your map axes, and display your map.

figure
axesm globe
axis off
geoshow(A,R)
title('Blue Marble: Next Generation')

The layer used in this example is courtesy of NASA/JPL-Caltech.

Search the WMS database for layers that contain foundation.gtopo30 in the LayerName field. GTOPO30, a digital elevation model developed by the United States Geological Survey (USGS), has a horizontal grid spacing of 30 arc seconds.

gtopo30Layer = wmsfind('foundation.gtopo30');

Define a background color, specifying red, green, and blue levels.

oceanColor = [0 170 255];

Use the BackgroundColor and CellSize parameters of the wmsread function to set the background color and cell size of your retrieved map.

cellSize = 0.5;
[A,R] = wmsread(gtopo30Layer,'BackgroundColor',oceanColor, ...
    'CellSize', cellSize);

Open a figure window and set up a world map axes. Display your map with a title.

figure
worldmap world
geoshow(A,R)
title({'GTOPO30 Elevation Model',gtopo30Layer.LayerTitle})

Search for layers from the NASA Goddard Space Flight SVS Image Server. This server contains layers that have legend images. You can tell that legend images are available because the layers have content in the LegendURL field.

layers = wmsfind('svs.gsfc.nasa.gov','SearchField','serverurl');
serverURL = layers(1).ServerURL;
gsfc = wmsinfo(serverURL);

Find the layer containing urban temperature signatures and display the abstract:

urban_temperature = gsfc.Layer.refine('urban*temperature');
disp(urban_temperature.Abstract)

Big cities influence the environment around them. For example,
urban areas are typically warmer than their surroundings.
Cities are strikingly visible in computer models that simulate
the Earth's land surface. This visualization shows average 
surface temperature predicted by the Land Information System (LIS)
for a day in June 2001. Only part of the global computation
is shown, focusing on the highly urbanized northeast corridor
in the United States, including the cities of Boston, New York,
Philadelphia, Baltimore, and Washington.

Additional Credit:
NASA GSFC Land Information System (http://lis.gsfc.nasa.gov/)

Read and display the layer. The map appears with different colors in different regions, but without a legend it is not clear what these colors represent.

[A,R] = wmsread(urban_temperature);
figure
usamap(A,R)
geoshow(A,R)
title('Urban Temperature Signatures')
axis off

Investigate the Details field of the urban_temperature layer. This layer has only one structure in the Style field. The Style field determines how the server renders the layer.

urban_temperature.Details

ans = 

  struct with fields:

     MetadataURL: 'http://svs.gsfc.nasa.gov/vis/a000000/a003100/a003152/a003152.fgdc'
      Attributes: [1x1 struct]
     BoundingBox: [1x1 struct]
       Dimension: [1x1 struct]
    ImageFormats: {'image/png'}
     ScaleLimits: [1x1 struct]
           Style: [1x1 struct]
         Version: '1.3.0'

Display the Style field in the Command Window:

urban_temperature.Details.Style

ans = 

        Title: 'Opaque'
         Name: 'opaque'
     Abstract: [1x319 char]
    LegendURL: [1x1 struct]

Each Style element has only one LegendURL. Investigate the LegendURL:

urban_temperature.Details.Style.LegendURL

ans = 

    OnlineResource: [1x65 char]
            Format: 'image/png'
            Height: 90
             Width: 320

Download the legend URL:

url = urban_temperature.Details.Style.LegendURL.OnlineResource

The URL appears in the command window:

url =

http://svs.gsfc.nasa.gov/vis/a000000/a003100/a003152/temp_bar.png

Display the legend image using the image command and set properties, such that the image displays with one-to-one, screen-to-pixel resolution. This legend is simply an image of a colorbar, not the legend in MATLAB^® graphics.

temperatureLegend = webread(url);
figure('Color','white')
axis off image
set(gca,'units','pixels','position',...
[0 0 size(temperatureLegend,2) size(temperatureLegend,1)]);
pos = get(gcf,'position');
set(gcf,'position',...
[pos(1) pos(2) size(temperatureLegend,2) size(temperatureLegend,1)]);
image(temperatureLegend)

Now the map makes more sense. The regions toward the red end of the spectrum are warmer.

Steps 7–10 demonstrate how to capture the output from a map frame and append the legend. By appending the legend in this fashion, you avoid warping text in the legend image. (Legend text warps if you display the image with geoshow.)

First set your latitude and longitude limits to match the limits of your map and read in a shapefile with world city data:

latlim = R.LatitudeLimits;
lonlim = R.LongitudeLimits;
S = shaperead('worldcities','UseGeoCoords',true,...
    'BoundingBox',[lonlim(1) latlim(1);lonlim(2) latlim(2)]);

Determine the position of the current figure window. Vary the pos(1) and pos(2) 'Position' parameters as necessary based on the resolution of your screen.

colValue = [1 1 1];
dimension = size(A,1)/2;
figure
set(gcf,'Color',[1,1,1])
pos = get(gcf, 'Position');
set(gcf, 'Position', [pos(1) pos(2) dimension dimension])

Display the map and add city markers, state boundaries, meridian and parallel labels, a title, and a North arrow:

usamap(A,R)
geoshow(A,R)
geoshow(S, 'MarkerEdgeColor', colValue, 'Color', colValue)

geoshow('usastatehi.shp', 'FaceColor', 'none',...
   'EdgeColor','black')
mlabel('FontWeight','bold')
plabel('FontWeight','bold')
axis off
title('Urban Temperature Signatures', 'FontWeight', 'bold')

for k=1:numel(S)
   textm(S(k).Lat, S(k).Lon, S(k).Name, 'Color', colValue,...
       'FontWeight','bold')
end

lat = 36.249;
lon = -71.173;
northarrow('Facecolor', colValue, 'EdgeColor', colValue,...
    'Latitude', lat, 'Longitude', lon);

Display the map and legend as a single, combined image:

f = getframe(gcf);
legendImg = uint8(255*ones(size(temperatureLegend,1),size(f.cdata,2),3));
offset = dimension/2;
halfSize = size(temperatureLegend, 2)/2;
legendImg(:,offset-halfSize:offset+halfSize-1,:) = temperatureLegend;
combined = [f.cdata; legendImg];
figure
pos = get(gcf,'position');
set(gcf,'position',[10 100 size(combined,2) size(combined,1)])
set(gca,'units','normalized','position', ...
   [0 0 1 1]);
image(combined)
axis off image

Another way to display the map and legend together is to burn the legend into the map at a specified location. To view the image, use the image command, setting the position parameters such that there is a one-to-one pixel-to-screen resolution. (Legend text warps if the image is displayed with geoshow.)

A_legend = A;
A_legend(end-size(temperatureLegend,1):end-1,...
   end-size(temperatureLegend,2):end-1,:) = temperatureLegend;
figure
image(A_legend)
axis off image
set(gca,'Units','normalized','position',...
   [0 0 1 1]);
set(gcf,'Position',[10 100 size(A_legend,2) size(A_legend,1)]);
title('Urban Temperature Signatures', 'FontWeight', 'bold')

Combine the map and legend in one file, and then publish it to the Web. First write the images to a file:

mkdir('html')
imwrite(A_legend, 'html/wms_legend.png')
imwrite(combined, 'html/combined.png')

Open the MATLAB Editor, and paste in this code:

%%
% <<wms_legend.png>>

%%
% <<combined.png>>

Add any other text you want to include in your published document. Then select one of the cells and choose File > Save File and Publish from the menu.

Obtain geographic coordinates that are coincidental with the image in the file boston.tif.

filename = 'boston.tif';
info = georasterinfo(filename);
R = info.RasterReference;
proj = R.ProjectedCRS;
[latlim,lonlim] = projinv(proj, ...
                          R.XWorldLimits,R.YWorldLimits);

Convert the geographic limits to Web Mercator. To do this, first create a projcrs object by specifying the EPSG:3857 coordinate system. Then, project the coordinates using projfwd. To obtain the imagery in this coordinate reference system, you need to use WMSMapRequest and WebMapServer since wmsread only requests data in the EPSG:4326 system.

p = projcrs(3857);
[x,y] = projfwd(p,latlim,lonlim);
xlimits = [min(x) max(x)];
ylimits = [min(y) max(y)];

Calculate image height and width values for a sample size of 5 meters.

metersPerSample = 5;
imageHeight = round(diff(ylimits)/metersPerSample);
imageWidth  = round(diff(xlimits)/metersPerSample);

Re-compute the new limits.

yLim = [ylimits(1), ylimits(1) + imageHeight*metersPerSample];
xLim = [xlimits(1), xlimits(1) + imageWidth*metersPerSample];

Find the USGS National Map from the WMS database and select the Digital Ortho-Quadrangle layer.

doqLayer = wmsfind('usgsnaipplus','SearchField','serverurl');
doqLayer = wmsupdate(refine(doqLayer,'image'));
doqLayer = refine(doqLayer,'USGSNaip','SearchField','abstract');

Create WebMapServer and WMSMapRequest objects.

server  = WebMapServer(doqLayer.ServerURL);
numberOfAttempts = 50;
attempt = 0;
request = [];
while(isempty(request))
    try
        request = WMSMapRequest(doqLayer, server);
    catch e 
        attempt = attempt + 1;
        fprintf('%d\n', attempt)
        if attempt > numberOfAttempts
            throw(e)
        end
    end
end

Use WMSMapRequest properties to modify different aspects of your map request, such as map limits, image size, and coordinate reference system code. Set the map limits to cover the same region as found in the boston.tif file.

request.CoordRefSysCode = 'EPSG:3857';
request.ImageHeight = imageHeight;
request.ImageWidth  = imageWidth;
request.XLim = xLim;
request.YLim = yLim;

Request a map of the ortho-imagery in Web Mercator coordinates.

A_PCS = getMap(server, request.RequestURL);
R_PCS = request.RasterReference;

Obtain a map for the same region, but in EPSG:4326 coordinates.

request.CoordRefSysCode = 'EPSG:4326';
request.Latlim = latlim;
request.Lonlim = lonlim;
A_Geo = getMap(server, request.RequestURL);
R_Geo = request.RasterReference;

Read in Boston place names from a shapefile and overlay them on top of the maps. Convert the coordinates of the features to Web Mercator and geographic coordinates. The point coordinates in the shapefile are in meters and Massachusetts State Plane coordinates, but the GeoTIFF projection is defined in survey feet.

S = shaperead('boston_placenames');
x = [S.X]*unitsratio('sf','meter');
y = [S.Y]*unitsratio('sf','meter');
names = {S.NAME};
[lat, lon] = projinv(proj, x, y);
[xPCS, yPCS] = projfwd(p, lat, lon);

Project and display the ortho-imagery obtained in EPSG:4326 coordinates using geoshow.

figure
axesm('mercator','MapLatLimit',latlim,'MapLonLimit',lonlim)
geoshow(A_Geo,R_Geo)
textm(lat, lon, names, 'Color',[0 0 0], ...
    'BackgroundColor',[0.9 0.9 0],'FontSize',6);
axis tight
title({'USGS Digital Ortho-Quadrangle - Boston', ...
    'Geographic Layer'})

Display the ortho-imagery obtained in Web Mercator coordinates.

figure
mapshow(A_PCS, R_PCS);
text(xPCS, yPCS, names, 'Color',[0 0 0], ...
    'BackgroundColor',[0.9 0.9 0],'FontSize',6,'Clipping','on');
axis tight
title({'USGS Digital Ortho-Quadrangle - Boston', 'Web Mercator Layer'})

Display the image from boston.tif for comparison.

figure
mapshow(filename)
text(x, y, names, 'Color',[0 0 0], ...
    'BackgroundColor',[0.9 0.9 0],'FontSize',6,'Clipping','on');
axis tight
title({filename, proj.Name})