Retrieving Your Map

Ways to Retrieve Your Map

To retrieve a map from a WMS server, use the function wmsread or, in a few specific situations, the WebMapServer.getMap method. Use the getMap method when:

  • Working with non-EPSG:4326 reference systems

  • Creating an animation of a specific geographic area over time

  • Retrieving multiple layers from a WMS server

In most cases, use wmsread to retrieve your map. To use wmsread, specify either a WMSLayer object or a map request URL. Obtain a WMSLayer object by using wmsfind to search the WMS Database. Obtain a map request URL from:

  • The output of wmsread

  • The RequestURL property of a WMSMapRequest object

  • An Internet search

The map request URL string is composed of a WMS server URL with additional WMS parameters. The map request URL can be inserted into a browser to make a request to a server, which then returns a raster map.

Understanding Coordinate Reference System Codes

When using wmsread, request a map that uses the EPSG:4326 coordinate reference system. EPSG stands for European Petroleum Survey Group. This group, an organization of specialists working in the field of oil exploration, developed a database of coordinate reference systems. Coordinate reference systems identify position unambiguously. Coordinate reference system codes are numbers that stand for specific coordinate reference systems.

EPSG:4326 is based on the 1984 World Geodetic System (WGS84) datum and the latitude and longitude coordinate system, with angles in degrees and Greenwich as the central meridian. All servers in the WMS Database, and presumably all WMS servers in general, use the EPSG:4326 reference system. This system is a requirement of the OGC® WMS specification. If a layer does not use EPSG:4326, Mapping Toolbox™ software uses the next available coordinate reference system code. The Mapping Toolbox does not support automatic coordinate reference systems (systems in which the user chooses the center of projection). For more information about coordinate reference system codes, please see the Spatial Reference Web site.

Retrieving Your Map with wmsread

NASA's Blue Marble Next Generation layer shows the Earth's surface for each month of 2004 at high resolution (500 meters/pixel). Read and display the Blue Marble Next Generation layer.

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

    nasa = wmsfind('nasa', 'SearchField', 'serverurl');
    
  2. 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');
    
  3. 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 referencing matrix that ties A to the EPSG:4326 geographic coordinate system. The geographic limits of A span the full latitude and longitude extent of layer.

  4. 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')

Setting Optional Parameters

The wmsread function allows you to set many optional parameters, such as image height and width and background color. This example demonstrates how to view an elevation map in 0.5-degree resolution by changing the cell size, and how to display the ocean in light blue by setting the background color. For a complete list of parameters, see the wmsread reference page.

  1. The ORNL DAAC WMS server hosts a wide variety of layers. Search the WMS Database for layers that contain the string 'ornl.gov' in the ServerURL field.

    ornlLayers = wmsfind('ornl.gov', 'SearchField', 'serverurl');
    
  2. GTOPO30, a digital elevation model developed by the United States Geological Survey (USGS), has a horizontal grid spacing of 30 arc seconds. Refine your search to include only layers that contain the string 'gtopo30' in the LayerName and LayerTitle fields.

    gtopo30Layer = refine(ornlLayers, 'gtopo30');

    The refined search returns one layer.

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

    oceanColor = [0 170 255];
    
  4. 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);
    
  5. 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})
    

Adding a Legend to Your Map

A WMS server renders a layer as an image. Without a corresponding legend, interpreting the pixel colors can be difficult. Some WMS servers provide access to a legend image for a particular layer via a URL that appears in the layer's Details.Style.LegendURL field. (See the WMSLayer.Details reference page for more information.)

Although a legend provides valuable information to help interpret image pixel colors, only about 45% of the servers in the WMS database contain at least one layer with an available legend. Less than 10% of the layers in the WMS database contain a legend, but nearly 80% of the layers in the database are on the columbo.nrlssci.navy.mil server. This server always has empty LegendURL fields. You cannot use wmsfind to search only for layers with legends because the database does not store this level of detail. You must update a layer from the server before you can access the LegendURL field.

This example demonstrates how to create a map of surface temperature, and then obtain and display the associated legend image:

  1. 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);
    
  2. 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/)
  3. 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
    

  4. 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 = 
    
         MetadataURL: [1x65 char]
          Attributes: [1x1 struct]
         BoundingBox: [1x1 struct]
           Dimension: [1x1 struct]
        ImageFormats: {'image/png'}
         ScaleLimits: [1x1 struct]
               Style: [1x1 struct]
             Version: '1.1.1'

    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
  5. Download the legend URL and save it:

    url = urban_temperature.Details.Style.LegendURL.OnlineResource
    urlwrite(url, 'temp_bar.png');
    

    The URL appears in the command window:

    url =
    
    http://svs.gsfc.nasa.gov/vis/a000000/a003100/a003152/temp_bar.png
    
  6. Display the legend image using the image command and set properties, such that the image displays with one-to-one, screen-to-pixel resolution.

    legend = imread('temp_bar.png');
    figure('Color','white')
    axis off image
    set(gca,'units','pixels','position',...
       [0 0 size(legend,2) size(legend,1)]);
    pos = get(gcf,'position');
    set(gcf,'position',...
       [pos(1) pos(2) size(legend,2) size(legend,1)]);
    image(legend)
    

    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.

  7. 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,lonlim] = limitm(A,R);
    S = shaperead('worldcities', 'UseGeoCoords',true,...
        'BoundingBox',[lonlim(1) latlim(1);lonlim(2) latlim(2)]);
    
  8. 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])
    
  9. 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);
    

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

    f = getframe(gcf);
    legendImg = uint8(255*ones(size(legend,1), size(f.cdata,2), 3));
    offset = dimension/2;
    halfSize = size(legend, 2)/2;
    legendImg(:,offset-halfSize:offset+halfSize-1,:) = legend;
    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
    

  11. 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(legend,1):end-1,...
       end-size(legend,2):end-1,:) = legend;
    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')
    
  12. 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.

Retrieving Your Map with WebMapServer.getMap

The WebMapServer.getMap method allows you to retrieve maps in any properly defined EPSG coordinate reference system. If you want to retrieve a map in the EPSG:4326 reference system, you can use wmsread. If you want to retrieve a layer whose coordinates are not in the EPSG:4326 reference system, however, you must use the WMSMapRequest class to construct the request URL and the WebMapServer.getMap method to retrieve the map. This example demonstrates how to create maps in "Web Mercator" coordinates using the WMSMapRequest and WebMapServer classes. The Web Mercator coordinate system is commonly used by web applications.

The USGS National Map provides ortho-imagery and topography maps from various regions of the United States. The server provides the data in both EPSG:4326 and in Web Mercator coordinates, as defined by EPSG codes EPSG:102113, EPSG:3857. For more information about these codes, see the Spatial Reference Web site.

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

    filename = 'boston.tif';
    proj = geotiffinfo(filename);
    cornerLat = [proj.CornerCoords.Lat];
    cornerLon = [proj.CornerCoords.Lon];
    latlim = [min(cornerLat) max(cornerLat)];
    lonlim = [min(cornerLon) max(cornerLon)];
  2. Convert the geographic limits to Web Mercator. The EPSG:102113 coordinate system, is commonly known as "Web Mercator", or "Spherical Mercator projection coordinate system". This system uses the WGS84 semimajor axis length (in meters) but sets the semiminor axis length to 0. 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.

    mstruct = defaultm('mercator');
    mstruct.origin = [0 0 0];
    mstruct.maplatlimit = latlim;
    mstruct.maplonlimit = lonlim;
    wgs84 = wgs84Ellipsoid;
    mstruct.geoid = [wgs84.SemimajorAxis 0 ];
    mstruct = defaultm(mstruct);
    
    [x,y] = mfwdtran(mstruct, latlim, lonlim);
    xlimits = [min(x) max(x)];
    ylimits = [min(y) max(y)];
  3. 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);
    
  4. Re-compute the new limits.

    yLim = [ylimits(1), ylimits(1) + imageHeight*metersPerSample];
    xLim = [xlimits(1), xlimits(1) + imageWidth*metersPerSample];
  5. Find the USGS National Map from the WMS database and select the Digital Ortho-Quadrangle layer.

    nmap = wmsfind('nationalmap', 'SearchField', 'serverurl'); 
    doqLayer = refine(nmap, 'NAIP*ImageServer', 'SearchField', 'serverurl'); 
    doqLayer = doqLayer(1); 
    
  6. 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
    
  7. 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:102113';
    request.ImageHeight = imageHeight;
    request.ImageWidth  = imageWidth;
    request.XLim = xLim;
    request.YLim = yLim;
  8. Request a map of the ortho-imagery in Web Mercator coordinates.

    A_PCS = getMap(server, request.RequestURL);
    R_PCS = request.RasterRef;
  9. 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.RasterRef;
  10. 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] = mfwdtran(mstruct, lat, lon);
  11. Project and display the ortho-imagery obtained in EPSG:4326 coordinates using geoshow.

    mstruct.geoid = wgs84;
    mstruct = defaultm(mstruct);
    
    figure('Renderer','zbuffer')
    axesm(mstruct)
    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'})

  12. 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'})

  13. 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.PCS})
    

Was this topic helpful?